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Abstract 

Sterile neutrinos are SU{2) singlets that mix with active neutrinos via a mass matrix, its di- 
agonalization leads to mass eigenstates that couple via standard model vertices. We study the 
cosmological production of heavy neutrinos via standard model charged and neutral current ver¬ 
tices under a minimal set of assumptions: i) the mass basis contains a hierarehy of heavy neutrinos, 
ii) these have very small mixing angles with the active (flavor) neutrinos, iii) standard model parti¬ 
cles, including light (active-like) neutrinos are in thermal equilibrium. If kinematically allowed, the 
same weak interaction processes that produce active-like neutrinos also produce the heavier species. 
We introduce the quantum kinetic equations that describe their production, freeze out and decay 
and discuss the various processes that lead to their production in a wide range of temperatures 
assessing their feasibility as dark matter candidates. The final distribution function at freeze-out 
is a mixture of the result of the various production processes. We identify processes in which finite 
temperature colleetive exeitations may lead to the production of the heavy species. As a specific 
example, we consider the production of heavy neutrinos in the mass range < 140 MeV from 
pion decay shortly after the QCD crossover including finite temperature corrections to the pion 
form factors and mass. We consider the different decay channels that allow for the production of 
heavy neutrinos showing that their frozen distribution functions exhibit effects from “kinematic 
entanglement” and argue for their viability as mixed dark matter candidates. We discuss abun¬ 
dance, phase space density and stability constraints and argue that heavy neutrinos with lifetime 
r > 1/Hq freeze out of local thermal equilibrium, and eonjecture that those with lifetimes r <C 1/Hq 
may undergo cascade decay into lighter DM candidates and/or inject non-LTE neutrinos into the 
cosmic neutrino background. A comparison is made between production through pion decays with 
the production of non-resonant production via active-sterile mixing. 
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I. INTRODUCTION 


The present understanding of the large scale cosmological evolution is governed by un¬ 
known quantities, dark matter and dark energy [l|, where the commonly accepted wisdom 
points towards a dark matter candidate which is a weakly-interacting, cold thermal relic 
(WIMP) [^. Dark matter candidates alternative to the WIMP paradigm attract a fair 
amount of attention j^, due to several observational shortcomings of the standard cold 
dark matter cosmology (ACDM) at small scales. N-body simulations of cold dark mat¬ 
ter produce dark matter profiles that generically feature cusps yet observations suggest a 
smooth-core profile j^, [^(core-cusp problem). The same type of N-body simulations also 
predict a large number of dark matter dominated satellites surrounding a typical galaxy 
which is inconsistent with current observations (missing satellites problem). Both the 
missing satellites and core-cusp problem can be simultaneously resolved by allowing some 
fraction of the dark matter to be “warm” (WDM) jS-T^ with a massive “sterile” neutrino 
being one popular candidate 141-1 1^ - other examples include Kaluza-Klein excitations from 
string compactification or axions |3|, [lOj . The “hotness” or “coldness” of a dark matter can¬ 
didate is discussed in terms of its free streaming length. A/s, which is the cut-off scale in the 
linear power spectrum of density perturbations. Cold dark matter (CDM) features small 
(< pc) A/s that brings about cuspy prohles whereas WDM features A/s ~ few kpc which 
would lead to cored profiles. Recent WDM simulations including velocity dispersion suggest 
the formation of cores but do not yet reliably constrain the mass of the WDM candidate in 
a model independent manner since the distribution function of these candidates is also an 
important quantity which determines the velocity dispersion and thereby the free streaming 
length}^. 

For most treatments of sterile neutrino dark matter, a nonthermal distribution function 
is needed in order to evade cosmological bounds 2l|. The mechanism of sterile neutrino 

(see 


22 | 


production in the early universe through oscillations was originally studied in ref. 
also the reviews jl^) and in [l^ 23-27] sterile neutrinos are argued to be a viable warm 
dark matter candidate produced out of LTE in the absence (Dodelson-Widrow, DW) or 
presence (Shi-Fuller, SF) of a lepton asymmetry. Models in which a standard model Higgs 
scalar decays into a pair of sterile neutrinos at the electroweak scale (or higher) also yield an 
out-of-equilibrium distribution suitable for a sterile neutrino dark matter candidate j23-31|. 
Observations of the Andromeda galaxy with Chandra led to tight constraints on the (DW) 
model of sterile neutrinos j^. and more recently the observation of a 3.5 keV signal from the 
XMM Newton X-ray telescope has been argued to ^ due to a 7 keV sterile neutrino jH, 34 
a position which has not gone unchallenged j 


38 


The prospect of a keV sterile DM 

candidate continues to motivate theoretical and observational studies [13, 113, [111, 13-45 


In all of these production mechanisms, the assumption of a vanishing initial population is 
implemented and, as we discuss, there are a wide array of processes which are ignored with 
this simplihcation - a problem which is starting to become appreciated j46|. 

On the particle physics front, neutrino masses, mixing and oscillations are the hrst ev¬ 
idence yet of physics beyond the standard model. A signihcant world-wide experimental 
program has brought measurements of most of the parameters associated with neutrino 
mass 
future 


471 l48l with several significant remaining questions poised to be answered in the near 
' 49 I . Short baseline neutrino oscillation experiments (LSND, MiniBooNE) jsoj, 51 


present a picture of the neutrino sector which may require an additional sterile neutrino 


species of mass ~ leV [18|, l52|, l53|| but there remains tension with other experiments [54 
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and a definitive resolntion of these anomalies will reqnire farther experiments j55l-l60| . An 


interpretation of short baseline experiment anomalies as a signal for sterile nentrinos leads to 
a relatively light mass ~ eV for the sterile nentrino; however, it has been argued that sterile 
neutrinos with mass on the order of MeV or larger can decay and could also explain 
the short baseline anomalies or, alternatively, heavy sterile neutrinos produced through rare 

J. Well motivated proposals make the 


decay channels could also explain the anomaly 
case for a robust program to search for multiple heavy neutrinos 
of experiments including hadron colliders [6^68 


M in a wide range 


Furthermore, it has been argued that 
heavy sterile neutrinos in the mass range 100 — 500 MeV can decay nonthermally and enable 
evasion of cosmological and accelerator bounds j^. 

Several current experimental programs seek direct detection of sterile neutrinos : the 
KATRIN (Karlsruhe Tritium Neutrino Experiment) experiment 70,(3 searches for sterile 


neutrinos with masses up to < 18keV in tritium beta decay, the MARE experiment [3| (Mi¬ 
crocalorimeter Arrays for a Rhenium Experiment) explores the mass range < 2.5 keV in the 
beta decay of Rhenium 187, the ECHo experiment (Electron Capture Ho Experiment) (3| 
searches for sterile neutrinos in the mass range < 2keV in beta decay of Ho. Various re¬ 
cent proposals make the case for searches of heavy neutrinos at the Large Hadron Collider (3- 


68l | and current and future underground neutrino detectors may be able to probe dark matter 
candidates with ~ few MeV[3]- extensions beyond the standard model, the possibility of 
a hierarchy of heavy neutrinos would be natural: there is a (very wide) hierarchy of quark 
masses, charged lepton masses, and as is now clear a hierarchy of light massive neutrinos. 

Models of many dark matter components had been proposed recently[3]) these models 
provide the tantalizing possibility that the decay of one heavy component can seed the 
production of a lighter component. As we will discuss below, this possibility also arises if 
there is a hierarchy of heavy neutrinos. 


Motivation and Goals: As discussed above the study of sterile neutrinos with masses 
in a wide range, from few eV through keV up to several GeV is motivated from astrophysi- 
cal observations, cosmological simulations, terrestrial accelerators experiments and the new 
generation of experiments that will directly search for signals of sterile neutrinos. Most 
theoretical extensions beyond the Standard Model that provide mechanisms to generate 
neutrino masses invoke one or more generations of heavy “sterile” neutrinos. While the 
focus on sterile neutrinos as a dark matter candidate has been on the mass range of few 
keV (largely motivated by the cusp-vs-core and related problems of structure formation, and 
more recently by the possible detection of an X-ray line at ~ 7keV), if extensions beyond the 
standard model allow for a hierarehy of heavy neutrinos these may yield a mixture of warm, 
cold (and hot) dark matter. This possibility motivates us to explore possible mechanisms 
of production of heavier species of massive neutrinos and assess in particular cases their 
production, freeze-out and possible consequences of non-thermal distribution functions. 

Mixed dark matter described by several species of massive neutrinos with non-equilibrium 
distribution functions will certainly evade Lyman-a constraints 0 - 

Rather than proposing yet new models, the purpose of this work is twofold: 


• i): to understand the production and freeze out of heavy sterile-like neutrinos from 
standard model eharged and neutral eurrent interaetion vertices under a minimal set 
of assumptions, assessing their suitability as DM candidates. We argue that different 
processes with different kinematic channels are available in a wide range of temper- 
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atures and produce heavy neutrinos with different contributions to their distribution 
functions. The total distribution function is, therefore, a mixture of the contributions 
from the various different processes. 

• ii): to provide an explicit example within the context of the well understood effective 
held theory of weak interactions of pious including hnite temperature corrections. Pion 
decay is a production mechanism that is available shortly after the QCD crossover 
and produces heavy neutrinos from different kinematic channels. Therefore furnish an 
explicit example of the mixed nature of the distribution function. 

The minimal set of assumptions are the following: 

• We assume that the massless havor neutrinos helds of the standard model are related 
to the helds that create and annihilate the mass eigenstates as 

^ ^ T ^ ^ (^) (^■^) 

m=l,2,3 h=4,5--- 

where m = 1,2,3 refer to the light mass eigenstates that explain atmospheric and 
solar oscillation observations, h = 4, 5 ■ ■ ■ refer to heavy mass eigenstates. We assume 
a hierarchy of masses with Mm -C M^. The heavier mass eigenstates correspond to the 
various putative massive neutrinos with — leV to explain the LSND/MiniBooNE 
anomalies, or the — 50 — 80 MeV proposed to explain these anomalies via radiative 
decay [^. or M/j ~ 7 keV that could explain the X-ray line, or even — few MeV — 
GeV that could be a CDM candidate. We do not endorse nor adopt any particular 
extension beyond the Standard Model with particular mass generation mechanisms. 
We only assume the relationship (II.ip and the existence of a hierarchy of very massive 
neutrinos with mass scales well separated from the atmospheric and solar ones. 

• We assume that \Uam\ — 0{1) and \Hah\ ^ \Uam\ although this is a feature of generic 
see-saw type mechanisms, we do do not invoke a particular mechanism. We then write 
the charged and neutral current weak interaction vertices in the mass basis and keep 
only the linear terms in Hah- We only consider the production of the heavy species 
with Mh from these weak interaction vertices in the mass basis. In this manner the 
production of heavy neutrinos is similar to the production of standard model neutrinos 
via charged and neutral current vertices if the kinematics allows for the particular 
process to produce a heavy neutrino mass eigenstate. 

• Our last assumption is that the light mass eigenstates (active-like) ^ 1 ^ 2,3 along with all 
the other standard model particles are in thermal equilibrium. We will assume that 
the relevant processes occurred during the radiation dominated era with T > 0.1 MeV 
at which standard model active neutrinos decouple. 


This minimal set of assumptions allows us to implement well understood standard meth¬ 
ods from quantum kinetic theory to describe the production of heavy neutrinos from weak 
interaction vertices and asses their distribution function when the production freezes and 
discuss hnite temperature corrections to the production processes. 

Here we do not consider the possibility that the masses of sterile neutrinos is a conse¬ 
quence of Yukawa coupling to the Standard Model Higgs, as this entails a particular model. 
Production of heavy neutrinos from scalar decay has been discussed in refs. [2^31. 7^ . 
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As a definite example we implement this program in the particular case of production of 
heavy neutrinos from pion decay after the QCD transition (crossover) within the well under¬ 
stood effective held theory of charged pion decay including hnite temperature corrections to 
the pion decay constant and mass. This production channel is one of the most ubiquitous 
sources of neutrinos in accelerator experiments. Furthermore, pion decay shortly after the 
hadronization transition (at t ~ 10 fis) in the Early Universe certainly leads to neutrino 
production just as in accelerator experiments. 

It is shown that this production mechanism yields heavy neutrinos from different chan¬ 
nels (/i, e) which freeze out with highly non-thermal distributions, and furnishes a dehnite 
example of “kinematic entanglement” with important cosmological consequences on their 
clustering properties. 

Brief Summary of Results: 

• In the mass basis the same standard model processes that lead to the production of 
active-like (light) neutrinos yield heavy (“sterile-like”) neutrinos if kinematically al¬ 
lowed. We identify several processes that lead to the cosmological production of heavy 
neutrinos during the radiation dominated era in a wide range of temperatures includ¬ 
ing 'jUm —)■ i^h whose inverse process describes the X-ray telltale of sterile neutrinos 
as well as the possibility of production from collective excitations in the medium. To 
leading order in the mixing matrix elements Hah we obtain the quantum kinetic equa¬ 
tions and give their exact solution in terms of gain and loss rates that obey detailed 
balance and can be calculated with standard model rules. We analyze the possibility 
of thermalization and argue that heavy neutrinos with lifetimes larger than I/Rq will 
freeze-out with non-equilibrium distribution functions. We establish the bounds from 
abundance, coarsed grained phase space density (Tremaine-Gunn) and stability for 
suitable DM candidates and discuss clustering properties of heavy neutrino species in 
terms of the distribution function after freeze-out. 

• We generalize the concept of mixed DM to encompass not only different species of 
heavy neutrinos, but also the different distributions functions of a single species of 
mass Mh from different production channels. We argue that heavy neutrinos pro¬ 
duced via standard model charged and neutral current interactions are kinematically 
entangled with leptons produced in the same reaction. This kinematic entanglement 
leads to different distribution functions for different channels of production of a given 
species, some colder than others depending on the mass of the lepton, quantihed by the 
equation of state in each case. The total distribution function is therefore a mixture of 
all the different distribution functions from each channel. If the heavy neutrino does 
not thermalize (and those with lifetimes larger than 1/Hq do not), its distribution 
function at freeze-out exhibit memory of this kinematic entanglement. This memory 
is manifest in the equation of state, velocity dispersion, coarse grained phase space 
density and free streaming length. 

• As a specihc example, we study the production of heavy neutrinos after the QCD 

transition (crossover) from pion decay within the effective held theory description 
of weak interactions of pions including hnite temperature corrections to the pion 
decay constant and mass in a large kinematically allowed window Mh ^ 140 MeV. 
Specihcally we study in detail the distribution function from both available channels, 
TT —)■ ; TT —)■ evh, which furnishes a clear example of the “kinematic entanglement” 
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with the corresponding charged lepton. The distribution function from the /i channel 
is distinctly colder than that from the e channel and the total distribution function in 
this case is a sum (mixture) of both components, with concentrations depending on 
the ratio of mixing matrix elements \Hah\‘^ for each channel. We assess heavy neutri¬ 
nos produced via this mechanism as possible DM candidates by analyzing the allowed 
region of parameters that fulhlls the abundance, phase space density and stability 
constraints, highlighting how the phenomenon of kinematic entanglement is manifest 
along the boundaries of the allowed regions. We also study the equation of state (ve¬ 
locity dispersion) and free streaming length (cutoff in the power spectrum) and show 
that both interpolate between the colder (/r-channel) and warmer (e-channel) compo¬ 
nents as a function of Mh and the ratio showing explicitly the mixed 

nature of the distribution function. 

• We conjecture that heavy neutrinos with lifetimes shorter than 1 / Hq produced during 
the radiation dominated era may decay in a cascade into active-like neutrinos well 
after Big Bang Nucleosynthesis, providing a late injection of neutrinos out of thermal 
equilibrium into the cosmic neutrino background, and possibly, into lighter (but still 
heavy) and stable(r) neutrinos that could be DM candidates after matter radiation 
equality. 


II. PRODUCTION AND FREEZE-OUT: QUANTUM KINETICS 

Sterile neutrinos are SU{2) weak singlets that do not interact via standard model weak 
interactions, they only couple to the massless, flavor active neutrinos via an off diagonal mass 
matrix. This is the general description of sterile neutrinos, different models propose different 
types of (see-saw like) mass matrices. For our purposes the only relevant aspect is that the 
diagonalization of the mass matrix leads to mass eigenstates, and these are ultimately the 
relevant degrees of freedom that describe dark matter in its various possible realizations, cold, 
warm or hot. Cold and warm dark matter would be described by heavier species, whereas the 
usual (light) active-like neutrinos could be a hot dark matter component depending on the 
absolute value of their masses. As mentioned above in this article we focus on the production 
of heavy neutrinos from standard model couplings: in the mass basis, heavy neutrinos couple 
via standard model charged and neutral current interactions albeit with very small mixing 
matrix elements. If kinematically allowed, the same processes that produce active neutrinos 
will produce heavy (sterile-like) neutrinos with much smaller branching ratios. 

We consider the standard model augmented by neutrino masses under the assumption of 
a hierarchy of heavy neutrinos with masses much larger than those associated with the (light- 
active like) atmospheric and solar neutrinos. Upon diagonalization of the mass matrix, the 
active (massless) flavor neutrino helds are related to the neutrino helds that create/annihilate 
mass eigenstates as 


^o,l(^) ^ ^ T ^ ^ (II.1) 

m=l,2,3 h=4,5--- 

where a = e,/i, r, m = 1, 2, 3 refer to the light (atmospheric, solar) mass eigenstates with 
massess Mm and h = 4, 5 ■ • • refer to heavy mass eigenstates of mass Mh- The form of 
the mixing matrix is taken in a way that is model independent in that no assumptions are 
made as to the particular source of the neutrino masses and mixing. No specihcation of 
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the number of Dirac or Majorana mass terms nor their source is made. In particular, we 
do not assume that any Dirac mass terms are generated through Yukawa couplings to the 
Standard Model Higgs as this constitutes a choice of one particular model. We focus solely 
on production processes from standard model charged and neutral interactions. 

As discussed above our main assumptions are 

\Uam\ ^ Oil) > \H^h\ ; (II.2) 

and we assume that charged and neutral vector bosons, quarks, leptons and light 1 ^ 1 , 2,3 neu¬ 
trinos are all in thermal equilibrium. This assumption is justified with the usual arguments 
that the standard model interaction rates are much larger than the expansion rates down 
to T ~ 1 MeV when active-like neutrinos decouple from the plasma. To present the argu¬ 
ments in the simplest terms, we will consider vanishing chemical potentials for all relevant 
species, a posteriori, it is straightforward to include lepton asymmetries with chemical po¬ 
tentials. Furthermore, we will not distinguish between Dirac and Majorana neutrinos as the 
difference typically results in factors 2 in various transition probabilities. Lastly, we will 
not consider the possibility of CP-violation, implying that forward and backward (gain and 
loss) transition probabilities are the same. All of these assumptions can be relaxed with the 
concomitant complications which will be addressed elsewhere. 

In the mass basis the full Lagrangian density is 

C- = C,mSM + 5 (II-3) 

fc=4,5--- 

where CmSM is the Standard Model Lagrangian augmented with a diagonal neutrino mass 
matrix for the active like neutrinos Um of masses (m = 1, 2, 3), in this Lagragian density 
the charged and neutral interaction vertices are written in terms of the neutrino mass eigen¬ 
states with the mixing matrix Uam, ^oh is the free held Lagrangian density for the heavier 
neutrinos Uh of masses Mh (h = 4, 5 ■ ■ ■) and to linear order in Hah 

c, = -Any 5^ 5^ + h.c 

' a=e,iJ,,T h=4,5--- 

- Y E + h.c., (II.4) 

where ^ 

Q:=e,/i,r 

with Vl = (1 ~ 75)/2- From now onwards, we refer to heavy neutrinos instead of “sterile” 
neutrinos, because unlike the concept of a sterile neutrino, heavy neutrinos do interact with 
standard model degrees of freedom via charged and neutral current vertices. The new mass 
eigenstates will undoubtedly contribute to the tightly constrained “invisible width” of the 
ji^; however, the contribution from the heavier neutrinos is suppressed with respect to 
the light active-like neutrinos by small branching ratios Br ~ \Hih\‘^/\Uim\^ -C 1 and the 
tight constraints on the number of neutrinos contributing to the width of the may be 
evaded by very small matrix elements, which is, in fact, one of the underlying assumptions 
of this work. 
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The strategy is to pass to the Hamiltonian 


H = HmSM + H^{yh) + Hj = + Hi (II.6) 

where is the total Hamiltonian of the Standard Model in the mass basis of the light 
(active-like) neutrinos plus the free field Hamiltonian of the heavy neutrinos h = 4, 5 • • • and 
Hi is the interaction Hamiltonian obtained from £/ in (111.4^ . We now pass to the interaction 
picture of H^^'> wherein 

//,(()= (II.7) 

namely we will obtain the quantum kinetic equations to leading order in gw\Hah\ ^ Qw but 
in principle to all orders in the weak interaction coupling (more on this issue below). 
In terms of the interaction vertices flll.4p we can obtain the quantum kinetic equations for 
production of massive neutrinos to order {\Hah\Y from standard master (gain-loss) equations 
of the form 

dnh{q]t) _ dnh{q-,t ) , dnh{q-,t ) , 

dt ~ dt dt ^ ’ 

where nh{q]t) is the distribution function of the heavy neutrino. This is calculated for each 
process: the gain term describes the increase in the population nh{q]t) by the creation of a 
heavy state and the loss by the annihilation, the best way to understand the calculational 
aspects is with a few examples. 


A. Setting the stage in Minkowski space-time. 


We begin by describing the main aspects in Minkowski space time to highlight im¬ 
portant concepts, and generalize the formulation to an expanding cosmology in the next 
subsection. 

To begin with consider a temperature ^ T -C Tew where Tew — IbOGeV is the 

temperature of the electroweak transition. In this temperature range the massive vector 
bosons (described as three physical degrees of freedom in unitary gauge) are populated in 
thermal equilibrium in the plasma with the Bose-Einstein distribution functions. If the mass 
of the heavy neutrino is + rua < Mw,z fhe massive vector Bosons can decay into the 
massive neutrino, thereby contributing to the gain term. For example, take —)■ l^T^h] 

each la and each Uh constitute different decay channels which lead to a gain contribution 
whereas the loss contribution is the inverse process l^^h W^. The gain and loss terms 
are calculated by obtaining the total transition probability per unit time to a particular 
channel where an explicit calculation is detailed in appendix (j^. For h is 

straightforward to hnd 


271 


d^k\Mfi\^ 


dnh{q;t) I 

dt 2,Eh{q) j {271^2Ew{k)2Ea{j)) 

X 5{Ew{k) - Ea{p) - Eu{q)) 

dnh{q;t ), 


nw{k){l - n„(p))(l - nh{q-,t)) 


(11.9) 


277 


d^k \Mfi\ 


dt 2Eh{q) J {271)^2 Ew{k)2Ea{p) 

X 5{Ew{k) - Ea{p) - Eh{q)), 


{1 + nw{k))na{p)nh{q;t) 


( 11 . 10 ) 
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where p = \k — ^ and 


( 11 . 11 ) 


nw{k) 


1 

(,Ew{k)/T _ ^ 


Uaip) 


1 

(}Ea{p)/T _j_ ’ 


and \Aifi\‘^ is the nsnal transition matrix element for W~^ —)■ snmmed over the three 

polarizations of the vector Boson, snmmed over the helicity of the charge lepton and snmmed 
over the helicity of the heavy neutrino; obviously |A^/iP oc Therefore the quantum 

kinetic equation flll.SD becomes of the form 

= r^(g)(l - nh{q; t)) - T>{q)nh{q; t) (11.12) 

where the gain and loss rates are 


r<(g) 


27r 

‘iEhiq) 




Ott r \ M 

^ J (2.nEl(mEM (IM4) 

Because the W, are in thermal equilibrium the gain and loss rates obey the detailed balance 
condition 

r<(g)e-®^'*('?)/^ = r>(g), ( 11 . 15 ) 

which can be conhrmed straightforwardly from the explicit expressions fill. 13IIII. 141) using 
the energy conserving delta functions and the relations 

1 + nB{E) = e^/^nB{E) ; 1 - nF^E) = e^l^npiE) (11.16) 

where ub^f cU'e the Bose-Einstein and Fermi-Dirac distribution functions respectively. 

A similar exercise yields the quantum kinetic equation for the process —)■ and 

the inverse process, it is of the same form as flII.12p now with 


r<(g) 


27r 

2E/,(g) 


d^k\Mgi\^ 

{2n)^2Ez{k)2Em{p) 


nz{k){l-nm{p)) 6{Ez{k) - Em{p) - Eh{q)) 


(11.17) 


r>(g) 


2ti 

‘^Eh{q) 


d^k\Mfi\^ 

{2nf2Ez{k)2Em{p) 


{l+nz{k))n^{p)5{Ez{k)-E^{j))-Eh{q)) , 


(11.18) 


where now \A4fi\‘^ oc g‘^\Hhm\‘^ is the transition matrix element for the process —)• 

from the charged current vertex in fill.41) . Again because the vector Boson and the active-like 
neutrino Um are in thermal equilibrium, the gain and loss rates obey the detailed balance 
condition fill.151) . 

As highlighted above, the interaction picture corresponds to considering the usual stan¬ 
dard model vertices with light neutrinos Um to all orders in charged and neutral current 
weak interactions, as well as, in principle to all orders in strong interactions. This feature, 
in fact, leads to many different production channels of heavy neutrinos, provided the kine¬ 
matics allows the particular channels. In thermal equilibrium the typical energy of quarks 
or leptons whose masses are -C T is (E) ~ 3.15 T so that for Mh < T various higher order 
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processes are available; for example, f f ^ ^ i^m^h and its inverse process where /, / 

refer to quark-antiquark or lepton-antilepton or VmVm and similar charged current processes 
involving either quarks or charged leptons. If T ~ M\y the intermediate are 

on-shell, with nz{E),nw{E) ~ 1 for E T and the process is “resonantly” enhanced, the 
gain and loss terms are of the form 


dnh{q;t) I 

27r 


dt 

‘2Eh{q) ^ 

/ {27i)^2Ef{ki)2Ej{k2)2Em{p) 


X S{Ef{k,: 

I Ej{k2) - Em{p) - Eh{q)) 

dnh{q-,t) 1 

271 

r d^kid^k2\Mfi\^ 

dt 1'°"" 

2Eh{q) . 

1 {27i)^2Ef{ki)2Ej{k2)2Em{p) 


X S(E,(K 

) + Ej{k 2 ) - Em{p) - Eh{q)) 


nf{ki)nj{k2){l - n^(p))(l - nh{q]t)) 

(11.19) 

(1 - n/(fci))(l - nj{k2))nm{p)nh{q]t) 

( 11 . 20 ) 


where p = \ki + k 2 —^ and the contribution from the intermediate vector Boson is included in 
\A4fi\‘^, wherein the propagator must include the width (with the corresponding thermal 
contributions) because for T ~ Mz the intermediate state goes “on shell” and is enhanced. 
This results in a region in the phase space integrals with a large “resonant” contribution 
but with width ~ T^. A similar enhancement and treatment arises for processes of the form 
—)■ hT —)■ lai^h for T ~ My/ where the intermediate W becomes resonant. Obviously 
the difference (gain-loss) looks just like a typical Boltzmann equation, however, this is an 
equation for the production of the heavy neutrino Vh as the distribution functions of /, /, Vm 
are all in thermal equilibrium as per our previous assumption. For this contribution it fol¬ 
lows that oc g^\Hmh\‘^, for temperatures T <C Mz the Fermi limit implies that the 

production rate is oc G‘j;,T^Ei[Eh/T]\Hjnh\‘^ and similar contributions with charged current 
exchange, of the form G‘jpT^E 2 [Eh/T]\Hmh\‘^ on dimensional grounds, with Ei ^2 dimension¬ 
less functions of their arguments. A similar quantum kinetic equation describes the (gain) 
process / 1/2 ^ VF —>■ lai^h and its inverse (loss) process where /i ,/2 correspond to either 
up/down-type quarks or for charged current interactions with the concomitant change 
in oc It is clear that for these cases the general form of the quantum kinetic 

equation fill. 121) holds where again F<,F> obey the detailed balance condition flll.lSp . 

Another process that is important and relevant in cosmology is the production of via 
e'^e~h'rn ^h, the inverse process corresponds to the decay Uh e~^e~i'rn and a similar 
process mediated by neutral currents I'm^mEn —t Uh and the corresponding decay inverse 
process. In both cases the gain and loss terms are of the form 


dnh{q;t) 
dt 


I gain 


dnh{q; t) 
dt 


I loss 


27r 


d^kid^k 2 \Mfi 


2Eh{q) J {2TiY2Ei{ki)2E2{k2)2E^{p) 

X 5{E^{k^) + E2{k2) + E^{p) - Eniq)) 

2ti f d^kid^k2\Mlfi\‘^ 

2Eh{q) J {27ry2Ei{ki)2E2{k2)2E3{p) 

X nh{q; t) 5{Ei{ki) + .^ 2 (^ 2 ) + Es{p) - Eh{q)) , 


ni{ki)n2{k2)n^{p){l - nh{q;t)) 

( 11 . 21 ) 

(1 - ni(A;i))(l - n2{k2)){l - n^ip)) 


( 11 . 22 ) 
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where we read off the gain and loss rates 


r<(g) 


X 


27r F d^kid^k2\Aifi\‘^ 

2Eh{q) J {27r)^2Ei{ki)2E2{k2)2E3{p) 

6{E,{k,) + E2{k2) + E^{p) - Eh{q)) 


ni{ki)n2{k2)n^{p) 


r>(g) 


27r r 

2Eh{q) J {27r)^2Ei{ki)2E2{k2)2Es{p) ^ 

X S{E,{k,) + E2{k2) + Es{p) - Ef,{q)) , 


n2{k2)){l 


(11.23) 

nsip)) 

(11.24) 


where p = q — ki — k 2 and the labels 1, 2, 3 refer to the respective fermions either or 
Ujn for example (not to be confused with the labels for the light neutrinos Um)- We notice 
that as the temperature diminishes, setting the occupation factors uj = 0 in the loss term 
one recovers the decay rate of the heavy neutrino, this observation will be important in the 
discussion of thermalization and stability of the DM candidate below. 

It is clear from this discussion that in the mass basis, standard model charged and neu¬ 
tral current interaction vertices lead to production processes for heavy neutrinos that are 
similar to those of active-like neutrinos constrained by the usual kinematics. For example 
at temperatures T > 1 GeV, tan-lepton decay may lead to the production of heavy neu¬ 
trinos with masses up to ~ GeV. Heavy lepton decay is an available mechanism down to 
T ~ ~ lOOMeV and are processes that have not yet been studied in detail and clearly 

merit attention. 

Finally, at temperatures T > the production (gain) process —)• r'h is kinematically 

allowed, the inverse process Uh l^m is precisely the process conjectured to yield the X-ray 
line as a telltale of keV neutrinos. The corresponding Feynman diagrams for the gain and 
loss terms are depicted in fig. ([T]) For these processes we hnd 




7^m ^ 'jUmiloss) 

FIG. 1: yUm —)• r'h (gain) and the inverse process Vh 


T<{q) 


2 % 

2Eh{q) 


d^k\Msi\^ 

{27i)^2E^{k)2Em{p) 


n^{k)nm{p) 6{E^{k) Em{p) 


Ehid)) 


(11.25) 


F>(g) 


2ti 

‘2Eh{q) 


d^k\Mfi\^ 

{2'Kf2E^{k)2Em{,p) 


(1 -h n^{k)){l - nm{p))5{E^{k) E^{p) 


Eh{q)) , 

(11.26) 
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where the W — la loop is included in the \A4fi\‘^. How large can the production rate be?, at 
T -C Mw on dimensional grounds we expect (see discussion in section flllEp i 


^ ^ \HhiUim 

I 


(11.27) 


with a dimensionless function of the ratios Mh/T, Mi/T with a hnite limit for T S> 

For T ~ GeV this contribution to the production rate could be of the same order as that 
for non-resonant production (DW) at T ~ 150 MeVflsI. [ 2 ^, clearly motivating a deeper 
assessment of these processes. 

At T = 0 only the loss term flll.26p survives, and the corresponding decay rate has been 
obtained in ref. [78, 3|, this will be an important aspect discussed below. The gain process is 
actually kinematically allowed at temperatures T < because of the tail in the fermionic 


and blackbody distributions, although suppressed at lower temperatures. 

The important aspect of these latter processes is that whereas the gain rate T^{q) —)■ 0 
as the lepton and photon populations vanish at small temperature, the loss rate T^{q) does 
not vanish when the lepton populations vanish, the reason for this is that any population 
of Uh decays into the lighter leptons. These particular processes will be important in the 
discussion within the cosmological context. The radiative decay Uh —t i^ml is conjectured 
to be a telltale of the presence of “sterile” ( heavy) neutrinos. We then see that the inverse 
process produces the heavy (sterile-like) neutrinos at high enough temperature. 


B. Finite temperature corrections: 


There are important loop corrections at hnite temperature, self-energy corrections to the 
incoming and outgoing external “legs” as well as vertex corrections. There are also hnite 
temperature corrections to the mixing angles arising from self-energy loops, these tend to 


suppress the mixing matrix elements [23l. l80l . l81| therefore in medium the matrix elements 


Hfii —)■ Hhi^eff{T) and typically in absence of MSW resonances Hhi^eff{T) < Hm- An explicit 
example is given in section fllVp below (see eqn. flIV.2D b At high temperature there are hard 


thermal loop corrections to the self-energy of fermions and vector bosons [82 h 86[| that lead 
to novel collective excitations with masses oc gT where g is the gauge coupling. Photons 
and fermion-antifermion pairs form plasmon collective excitations with mass oc eT. For 
T > Tew the W, Z vector bosons do not acquire a mass via the Higgs mechanism because 
the ensemble average of the Higgs held vanishes, but they acquire thermal masses oc g^T 
akin to the plasmon collective excitations. Thermal masses for collective excitations of 
W, Z could open up kinematic windows for decay into above the electroweak transition. 
Plasmon collective excitations from photons in the medium can also produce via the 
electromagnetic process 7* —)■ via the Feynman diagrams displayed in hg. ([2]). These 

processes are similar to the mechanism of energy loss by plasmon decay into neutrinos in 
highly evolved massive stars 87H89| such as red giants and are also available in the Farly 
Universe. 

Although a priori these processes are subleading being suppressed by higher orders in 
the couplings, photons are populated all throughout the thermal history of the Universe, 
so these processes may contribute to production during a longer time scale as the leading 
processes described above. A similar possibility may be associated with plasmon collective 
excitations for above Tew resulting in the production of the heavier species at T > Tew, 


12 























this possibility merits deeper understanding, which is beyond the realm of this article. A 
follow up study will be reported elsewhere. 

The main point of this discussion is to highlight that in the mass basis, standard model 
interactions provide a wide variety of mechanisms to produce heavy neutrinos which are 
available at high temperatures in the early Universe. The final distribution function of a 
particular heavy neutrino species after freeze-out is a mixture of the distribution functions 
arising from the various production processes. A full assessment of a particular species as 
a DM candidate thus requires a thorough understanding of the various physical processes 
that lead to its production. 




FIG. 2: Gain processes: 7 * —>■ Vm^h 


C. Thermalization in Minkowski space-time 

An important aspect of the general form of the quantum kinetic equation flII.12j) is the 
linear dependence on the (time dependent) distribution function nh{q]t). This linearity is a 
consequence of keeping the lowest order term in \Hah\] \Ham\ and along with detailed balance 
has profound consequences: in Minkowski space time the quantum kinetic equation fill.121) 
leads to thermal equilibration of the heavy neutrino as the following argument shows. 

The solution of flII.12|) is 

= + ; 7(g) = r<(g) + r>(g), (11.28) 

Jo 

carrying out the time integral and using the detailed balance result fill.151) one finds 

nh{q; t) = neg{q) -f [nh{q; 0) - riegiq)] , (11.29) 

where neq{q) is the equilibrium (Fermi-Dirac) distribution function. Obviously 7 (g) is the 
rate of relaxation towards equilibrium, for t 3> l/ 7 (g) the distribution function is that of 
thermal equilibrium. Since 7 oc g‘^\H\‘^ one would be tempted to neglect the exponential 
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terms in flII.28D . however the exponent is secular, growing with time and becoming non- 
perturbative on a time scale t ~ I/ 7 . In Minkowski space time where the gain and loss rates 
are constant in time the distribution function will always reach thermal equilibrium at long 
time t 3> 1 / 7 . This is an important observation: detailed balance between the gain and loss 
term guarantees that at asymptotically long time t 3> I /7 the heavy neutrino thermalizes 
and its distribution function is figg. 

This is a general result in Minkowski space time which will be revisited below within the 
context of an expanding cosmology, where the relevant time scales are determined by the 
Hubble time scale. 


D. Production, freeze-out, LTE and decay with cosmological expansion: 

In an expanding cosmology described by an isotropic and homogeneous Friedmann- 
Robertson-Walker metric, and during a radiation dominated epoch the temperature varies 
with time T{t) oc l/a(t)[90[ where a{t) is the scale factor. The energy of a particle species 
of mass M measured locally by an observer is 

E(t) = + M-- ; p,[t) = ^ , (11.30) 

where pf(t),pc are the physical and comoving momenta respectively. The distribution func¬ 
tion of a particle species in local thermodynamic equilibrium (LTE) is (in absence of chemical 
potentials) 

n^{E{ty,t) = ^ ^ , (11.31) 

for Bose-Einstein and Fermi-Dirac respectively. The ratio 

Pc 

= y = 7^, 

to 

is constant during a radiation dominated era and Tq would be the temperature of the plasma 
today, related to the temperature of the cosmic radiation by accounting of the reheating of 
the photon gas whenever the number of the relativistic degrees of freedom changes j90[|. 
As the temperature drops during the expansion when M/T{t) 3> 1 the population of the 
massive species is strongly thermally suppressed (by setting vanishing chemical potential we 
assume that there is no conservation of this species). In LTE the calculation of the gain and 
loss rates carry over from Minkowski space-time by replacing the energies and momenta by 
E(t),pf{t) (IIL30p and the distribution functions of the species in LTE by flIL3ip . Although 
it is typical to separate the explicit time dependence in the rate equations from the time 
dependence of Pf{t), we will consider the rates and distribution functions as functions of 
Pc and the explicit time dependence of a{t). Thus the quantum kinetic equation in the 
cosmological setting becomes 

^ ~ nh{qc,t)) - r^{qc]t)nh{qc,t ), (11.33) 

where the gain and loss rates are obtained as in Minkowski space time replacing the momenta 
and energies by the local time dependent counterparts in the expanding cosmology and the 



(11.32) 
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(LTE) distribution functions for the various bosonic or fermionic species (with vanishing 
chemical potentials). 

As a consequence of the linearity of the quantum kinetic equation flll.33j) . which, in turn 
is a consequence of keeping only terms of 0{\Hah\‘^),0{\Hmh\^), the full quantum kinetic 
equation is a simple sum over all possible channels with total gain and loss rates 


all channels all channels 


(11.34) 


therefore the gain and loss rates in the quantum kinetic equation flll.33p are the total rates 

dirni) . 

Because for each channel the rates are calculated with distribution functions in 
(LTE) and obey the detailed balance condition this translates into 

t) = TlMc, t). (11.35) 


The general solution of flll.33p with the total gain and loss rates is 


nh{qc;t) = nh{qc;to) e 


- fto g- ftg 


(11.36) 


'to 


where the relaxation rate 


'y{qc, t) = t) + r>i(gc; t) 


neq{t) 


(11.37) 


with 

and we used the detailed balance result flll.35p . 

A hxed point solution with hh{qc]t) = 0 corresponds to freeze out. Although F^ and F> 
are related by detailed balance, the LTE distribution function rieq is not stationary, therefore 
is not a hxed point of the quantum kinetic equation^. However, the more relevant criterion 
for freeze-out is that the distribution function varies much more slowly than the expansion 
time scale, namely j^, 1 ^ . 

< Hit ). (11.39) 

nhiqc]t) 

As discussed above, the gain rate F”^ depends on the population of the particles whose 
decay or combination results in the production of the heavy neutrino, as a result the pro¬ 
duction rate diminishes during cosmological expansion eventually vanishing exponentially 
because of thermal suppression of the respective LTE distribution functions. Therefore if 


'jiqct') dt' -C 1 


'^0 


(11.40) 


1 


This is also the case for the Maxwell-Boltzmann distribution in a radiation dominated cosm 
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we can neglect the exponential terms in the solntion (III.36p . The remaining integral is 
hnite at asymptotically long time becanse the gain rate vanishes exponentially, leading to a 
non-equilibrium frozen distribution 


nh{qc]Oo) ~ nh{qc]to) + 


r£i(«c; «')*'■ 


(11.41) 


•to 


The freeze-ont condition is clearly achieved since flll.39p is fnlhlled with nh{qci t) = ^tbtidc] t) 
vanishing exponentially at snfficiently long time, whereas H{t) oc 1/t, and nh{qc] oo) 7 ^ 0 . 

When the condition flll.dOp is valid and the distribntion fnnction is linearly related to the 
prodnction rate as in flll.dip we can separate the contribntion from the different prodnction 
channels which will, generally, lead to different distribution functions, the total distribntion 
fnnction being a simple snm over all the different prodnction channels. Becanse in each 
prodnction channel the heavy nentrino is “kinematically entangled” with other leptons the 
distribntion fnnction at freeze-ont for each channel will reflect the different kinematics, which 
can be interpreted as a “memory” of the particnlar prodnction process. We will stndy a 
relevant example below. 

Althongh the gain and loss rates satisfy the detailed balance condition flll.35p the LTE 
distribntion fnnction is not a solntion of the qnantnm kinetic eqnation fill.331) becanse 
rieqit) 7 ^ 0. We can assess if and when a distribntion fnnction is very nearly in LTE by 
exploiting the relations flII.35IIII.37p to cast fill.331) in the form 


dnhjqclt) 

dt 


= nh{qc]t) -Ueqit) 


writing 


nhidc] t) = Ueqit) duhiqc, t ), 
and inserting into flll.42p one hnds 


(11.42) 

(11.43) 


Snh{qc;t) = 6nh{qc,to)e -e 


'*0 


with 


and 




(11.44) 

H(,t) 




H(t) = 1.664(r) ^ 2 X io=4(r)(M)%-i 


(11.45) 


(11.46) 


dnring radiation domination. As the temperatnre drops dnring expansion Ueq diminishes 
rapidly and we expect that the integral with neq{t) in flll.44p will be hnite at long time 
as the exponential snppression from neq{t) will overwhelm the pertnrbative growth of the 
integral of 7. Hence the condition that the distribntion fnnction becomes nearly LTE, 
becomes 


'y{qc,t')dt' > 1 . 


(11.47) 


'to 


at long time t < 1/Ho. If this condition is fnlhlled the full solution flIL36p must be considered 
as the integrals of 'y{qc, t) cannot be neglected. In this case the distribution function at long 
time may be described by (nearly) LTE, and the “memory” of the production process and the 
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kinematic “entanglement” characteristic of each prodnction channel is erased. Therefore, it 
is important to assess nnder what circumstances the condition (III.47P could be satished, since 
under such circumstances the heavy neutrino can thermalize with the rest of the standard 
model particles. 

As the temperature diminishes throughout the expansion history, the gain term is strongly 
suppressed by the thermal population factors because these processes entail the annihilation 
of particles in the initial state which are thermally populated (by assumption) (see the exam¬ 
ples described above). The loss rate will also vanish exponentially by thermal suppression at 
long time if it involves the annihilation of any thermal species, therefore for these processes 
the condition flll.dOp is expected to be fulhlled and the distribution function is expected to 
freeze out of LTE. 

However for processes in which the heavy neutrino decays into the hnal products the loss 
term does not vanish at low temperatures and the distribution function eventually decays, an 
example of this case is the loss term (111.241) describing the decay processes Uh —t e+e^z/^ or 
This has important implications: consider the quantum kinetic equation 
in the form flll.42p (note that 'j{qc]t) > 0) and an initial condition in the far past with 
nh{qc',to) = 0. At early times Uh grows h/i > 0 as the gain terms dominate and the heavy 
neutrinos are being produced. If at late times the population decays, namely < 0 this 
means that at some time the distribution function has reached LTE at which = 0, as 
the cosmological expansion continues the loss term dominates and the population decays. 
This simple analysis leads to the conclusion that if the lifetime of the heavy neutrino is 
smaller than the age of the Universe, meaning that it is now decaying, at some point in the 
past its distribution function reached LTE. Conversely, if the lifetime is much longer than 
the Hubble time 1/Hq then the distribution function has not reached LTE and the heavy 
neutrino is non-thermal. In other words, a heavy neutrino that is stable during the age of 
the Universe ~ 1/iLo and is a suitable DM candidate must be out of LTE. 

The analysis above indicates that when the lifetime of the heavy neutrino is much smaller 
than the age of the Universe, its distribution function must have passed through LTE during 
the evolution. The LTE solution neq{t) is not a true hxed point of the kinetic equation flH.33p 
because heq{Eh{f)) ^ 0 , however, we can ask what is the condition that ensures that the 
distribution function remains nearly in LTE if and when it reaches LTE. To answer this 
question we write 

nh{qc] t) = Uegit) (t) H- (11.48) 

where 6n^fl\t) <C neq{t)) etc, the expansion is in a small parameter to be quantihed a 


posteriori. Introducing this ansatz into flH.33j) we hnd 


neq{t) 

in terms of the relaxation time 


= (1 - neq{t)) 


r{qc] t) = 


Ml 


Hit) 


Tit)Ehif) 7(9cU) ’ 
1 


(11.49) 


(11.50) 


liqc^t) ’ 

we hnd that once LTE is attained, the distribution function will linger near LTE whenever 

Tit) 


Tiqc-,t)Hit) < 


M, 


(11.51) 


As the temperature drops, the contribution from T^ (gain) to 7 = r ^ is suppressed and the 
relaxation time becomes the lifetime of the heavy neutrino (loss via decay). This can be seen 
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for example, from the loss term for the process —)• with given by flII.24D with 

1,2,3 = m, m', m" and the discussion following it. Therefore the condition flll.Sip is actually 
equivalent to the statement of a lifetime much shorter than the Hubble time, consistently 
with the discussion above. If the heavy neutrino is a stable DM candidate its lifetime must 
be r > implying that its distribution function will not be in LTE. 


E. Stability and lifetime: 

To be a suitable DM candidate a heavy neutrino must feature a lifetime r > I/Hq, thus 
it remains to understand the decay channels for a hrmer assessment of the lifetime of heavy 
neutrinos, their suitability as DM candidates and the conditions for non-LTE distribution 
functions. 


The decay channels of a heavy neutrino were studied in refs. |57l. l58l. |79j, for example 
the purely leptonic and radiative channels: neutral current process —)■ or the 

charged current process Uh —t e~^e~h'rn and Uh t these are the inverse 

processes associated with the production processes —)■ i^h and 

^ r'/i respectively. At T = 0 the decay rates for z/^ —)■ e~^e~h'rn have 

been obtained in ref. jSTI . l58| , these are given by 


r(z//i e'^e Urn) — 3.5 x 10 ^\H, 


eh I 


Mh 

MeV 


K 


m" 

M 


where 


K [x] = (1 — (1 — 14a; — 2x^ — 12x^) + 24a;^(l — x^) In 


1 + (1 — 4a;) 

1 - (1 - 4a;)h2 


(11.52) 


(11.53) 


with iF —)■ 0 for —)■ 2me and iF —)■ 1 for 3> m-e. Therefore for this process the ratio 

of the lifetime to the Hubble time today 1/Hq is given by 

10"^^ f MeV\^ 

H,r{un ^ e+e-z.^) ^ \ w) ‘ 


\Heh?K 




The decay rate into active-like neutrinos mediated by neutral currents (not GIM 
fGlashow-Iliopoulos-Maiani) suppressed with sterile-like heavy neutrinos) is given by (see 


T{Uh = 3.5 X 10 ® ^ \H, 

therefore for this “invisible” channel we hud 


a.h 


Oi=e^fl^T 


Mh 

MeV 


HoT{Uh - 


10 


-13 


V 

Z-^a.= 


a=ey{i,T 




\Hr.h? V Mh 


(11.55) 


(11.56) 


The radiative decay Uh has been studied in ref. [58|, |78|, |79|, for heavy (sterile-like) 

neutrinos it is not GIM suppressed but is suppressed by one power of ttem with respect to 
the leptonic channels above and given by 


7 / -^h \ ^ ^ 

r(-'. ^ - io-4i^) 


ho/Ua.m I 


(11.57) 
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namely 


MeV 


5 


HoTi.l'h ll^m) ^ 


2 X 10 


-11 


Ea \HhaUa 


M. 


(11.58) 


For a Majorana neutrino there is an extra factor 2 multiplying the rates for the expressions 
above. For the heavy neutrino to be stable during the lifetime of the Universe and be a 
suitable DM candidate it must be that e+e^z/m) ; t ^ 

> 1 - 

Assuming that the \Hah\ are all of the same order and \Uam\ — 0{1) these estimates 
suggest that heavy neutrinos with 


Hho^WMh/MeVf < 10-13 


(11.59) 


feature lifetimes > 1/Hq and may be acceptable DM candidates. Interestingly, for Mh < 
1 MeV the “visible” leptonic decay channel Vh shuts off and the lifetime is 

dominated by the “invisible” channel Vh —t therefore evading potential bounds 

from the “visible” decays into Vl~ pairs. 

As per the discussion above the constrain that the heavy neutrino features a lifetime 
r > 1/Hq also implies that its distribution function is out of LTE. We will discuss possible 
interesting consequences of a heavy sterile neutrino with a lifetime much smaller than the 
Hubble time in a later section (see sectioi iV)) . 


F. Comparisons and caveats 


A formulation of the production rates of sterile neutrinos firmly based on the quantum 
field theory of neutrino mixing was introduced in ref. [1^. It relies on a see-saw type mass 
matrix with vanishing masses for the active neutrinos, large (Majorana) masses for sterile 
neutrinos on the diagonal and small off-diagonal matrix elements that mix the sterile and 
active neutrinos. This small (compared to the large Majorana masses) off-diagonal mixing 
sub-matrix is treated in a perturbative expansion. The authors in ref. 25|] obtain the quantum 
density matrix by tracing over the degrees of freedom of the standard model (assumed in 
thermal equilibrium) up to second order in the off-diagonal mixing, and, in principle, to all 
orders in the strong and weak interactions. The quantum master equation that describes 
the time evolution of the reduced density matrix is completely determined by correlation 
functions of the active neutrino including self-energy corrections, which is written in terms 
of spectral representations. The production rate of sterile neutrinos is obtained from the 
imaginary part of this self energy evaluated on the mass-shell of the sterile neutrinos, in 
principle to all orders in standard model couplings. The authors focus on temperature 
scales M\yz, and argue that sterile neutrinos with masses in the keV range are primarily 
produced in the temperature range T ~ 150 MeV. Furthermore, they only consider a “gain” 
term in this temperature regime neglecting the loss term. 


There are several differences between the approach of ref. [25j and the framework presented 
here: 1) we describe the production process with standard quantum kinetic equations for 
the mass eigenstates, without resorting to any particular model for the mixing mass ma¬ 
trix, however, similarly to our results are in principle valid to all orders in standard 
model couplings, but to second order in the (small) mixing matrix elements between the 
light (active-like) and heavy (“sterile-like”) neutrinos. 2) By going to the mass basis, we 
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recognize that standard model processes that produce active neutrinos via charged and neu¬ 
tral current vertices lead to the production of heavy neutrinos provided the kinematics is 
favorable and identify various processes available in a wide temperature region. 3) we con¬ 
sider both the gain and loss terms, giving an exact result for the non-equilibrium evolution 
of the distribution function, this is the result flll.36p . We discuss the important issue of 
thermalization vis a vis the constraints of stability of the DM candidate. 

The authors of ref. 2^ identihed the “one-loop” contributions to the active neutrino self 
energy, whose imaginary parts are precisely the contributions described by the quantum 
kinetic equations flII.13IIII.14IIII.17IIII.18p but neglected them by restricting their study to 
T -C Mw,z, whereas we argue that these contributions may be important (obviously a 
statement on their impact requires a detailed calculation, beyond the scope of this article). 
The production processes mediated by W,Z exchange in an intermediate state with gain and 
loss terms given by flII.191III.22P (and eqmvalent for charged current processes) correspond 
to the “two loops” contributions in ref. [2^, although we pointed out that at T ~ Mw^Mz 
these processes could be resonantly enhanced as the intermediate IT, Z propagators feature 
“on-shell” thermal contributions which are not suppressed at these temperatures. 4) We 
have identihed production processes from “collective excitations” in the medium, suggesting 
that hnite temperature corrections, leading to plasmon masses for the photon and for the 
W,Z bosons above Tew yield important production mechanisms. 5) we will argue below 
that different production channels with different kinematics yield different contributions to 
the distribution function. This is a result of “kinematic entanglement” between the heavy 
neutrino and the lepton produced in the reaction leading to distribution functions that may 
be colder for some channels and warmer for others depending on the mass of the lepton and 
the kinematics. In other words, the hnal distribution is a mixture of several components 
even for the same species. This point will be elaborated upon in more detail within a specihc 
example in the next sections. 

Caveats: As we discussed above at high temperature there are important self-energy 
corrections that must be assessed for a more reliable understanding of the gain and loss 
processes. Furthermore, our results are valid up to 0{\Hah\^) under the assumption that 
l-ffahP ^ I Damp, (and in principle to all orders in weak coupling), however this approxi¬ 
mation will break down if there are Mikheyev-Smirnov-Wolfenstein (MSW)j^, M, 92] res¬ 
onances in the medium, because near the resonance the effective mixing angle reaches vr/d. 
This is also a caveat in the approach of ref. because in this reference the quantum master 
equation has been obtained up to second order in the sterile-active mixing angle. 

Including the possibility of MSW resonances in the medium requires adopting a different 
framework that does not rely on \Hah\ ^ 1 but that would be hrmly based on quantum 
held theory out of equilibrium. Such approach could be based on effective held theory out 
of equilibrium as advocated in ref. We will report on this approach in a future study. 


III. COSMOLOGICAL CONSEQUENCES AND CONSTRAINTS. 

If there is a hierarchy of stable heavy neutrinos produced by the various mechanisms 
discussed above (we will comment later on unstable heavy neutrinos), each species with 
Mh 3> 1 eV will contribute as a non-relativistic DM component after matter radiation 
equality. Once the distribution function at freeze out is obtained, various cosmological 
consequences of the dark matter species can be assessed. In this section we gather together 
various quantities of cosmological relevance which are determined by basic properties of the 
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dark matter species: mass, number of intrinsic degrees of freedom, distribution function and 
freeze-out (or decoupling) temperature with the purpose of assessing particular species of 
heavy neutrinos as DM candidates once their distribution function is obtained. 

Let us define the total asymptotic “frozen” distribution function of a species of mass 
Mh as 

fh{qc) = nhiqc] oo), (III.l) 

where by f —)■ oo we mean a time sufficiently long that the distribution function satishes the 
freeze-out condition flll.39p . 

As discussed above after freeze-out the distribution function depends on the physical 
momentum qf{t) = qc/a{t) and temperature T{t) through the combination y = qf{t)/T{t) = 
qc/To where Tq is the temperature at which the species decoupled from the plasma (freeze- 
out) redshifted to today. Through the usual argument of entropy conservation it is related 
to the temperature of the cosmic microwave background today T^fl by 


-^ 7,0 ^9dh' 


(III.2) 


where Qdh is the number of ultrarelativistic degrees of freedom at the time when the particular 
species Uh decoupled (freezes). Then the number density energy density and partial 
pressure of species Uh is given by 








9u^ r i^/i^ 


9i^h f 




dy 


y'^fhiy) 


1 + 


+ M'j 

where is the number of internal degrees of freedom of the species and 

Xh{t) = Mh/T{t). 




(III.4) 


(III.5) 


(III.6) 


The contribution of each non-relativistic species Ph to today with Tq/M^ -C 1 is given 
by 


= 


^ - ( 7 f) (1.1.7) 

Using the relation flIII.2p and n^li^/pc = 1/25.67 eU we hnd 


g^^Mh ( Tq\^ 


2 ^ Mh / g^ 
61.7eU Jo 


dy y^fhiy) 
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Therefore, assuming that all of DM is in the form of various species of heavy neutrinos, it 
follows that 
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The fraction of ^dm contributed by the species Vh is given by 
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where we used = 0.1199[T7 . 

The equation of state for each species is 


wu,{Tm = 


puAt) 


rdy 


fh{y) 

m 2 
_|_ h 
T'^{t) 


/o°° dyy^sj: 


y2 I 

y r2(j) 


fh{y) 


(lll.ll) 


The equation of state yields information on how “cold” is the particular species and provides 
a benchmark to compare the equilibrium and non-equilibrium distribution functions. In 
particular yields a generalization of the effective adiabatic “speed of sound” for 

collisionless DM, namely, Vuf^it) = c\(T{t)) Pu^t) ; c\(T(t)) = Wy^(T(t)). 

In the non-relativistic limit S> T(t) ; x{t) —)■ cxo, the average primordial velocity 
dispersion for a species Vh is given by 




T{tf j dyy^fh{y) _ 3P, 
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which leads to the primordial velocity dispersion relation 


V = a 

' ’'h 


p^h 




_T{t) / Jdyy^fhiy) 


Mh V ^ I dyy'^fhiy) ’ 
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therefore for a collisionless non-relativistic species Uh the dispersion = Wui plays the role 
of the adiabatic sound speed. The equation of state (lIII.lip will be important in assessing 
the consequences of kinematic entanglement of distribution functions obtained by different 
production channels since it yields information on the “coldness” of the distribution. 


In analogy with a fluid description and following ref. [9^ we introduce the comoving free 
streaming wave vector for the species Vh similarly to the comoving Jeans wavevector, namely 




A-kG pm{t) 

(vm) 


aJ{t) 


(III.14) 


as discussed in ref. [9^ kfs{teq) describes the cutoff in the power spectrum for linear density 
perturbations, where teq is the time of matter-radiation equality. Since for non-relativistic 
particles 

<kt(o)> 


K-M = 


a^t) 


(111,16) 


where (K?^(0)) is the velocity dispersion today, it follows that 


f^fs,i'h{teq) \J^ijeq) 
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where 
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1/2 
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is the free streaming wavevector today. The free streaming scale that determines the length 
scale associated with the cutoff in the power spectrum is 




2tt 


kfs,Uh{^eq) 

with Vtuh'^ = + VtoMh^ = 0.14 we hnd^ 


— A/s,iv;i(0) (1 + Zeq) ^ 


(111.18) 


1 _ n / 2 y Ifdyy^fhiy) ^ mnnl 

A,.,„,(0) - 9,72 j y Jdvv-My) 

and (1 + ZeqY^'^ — 57. 

Up to constants of 0(1)) ^fs,uh{'^eq) is equivalent to the comoving distance traveled by 
a free streaming particle with average velocity (^)) from the time of matter-radiation 
equality until today jo^, as can be seen from the following argument. Consider a particle 

moving with velocity v{t) = \J the comoving distance traveled between teg and 
today at time t* ^ t^g is 


l{teq,t*) = 




a(t') 


dt' 


during a matter dominated cosmology when density perturbations grow 

a{t) = 

and using fllll.lSp for a non-relativistic species we hnd 

1 1/2 
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l{teq,t*) ~ 3 
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nuHS 


(1 + ZegY^"^ , 
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from which it follows that 

>^fs,i^Hiteq) ^ 1.7 l{teq,f) . (III.23) 

Phase-space density (Tremaine-Gunn) constraints are based on the result that the coarse 
grained phase space density may only decrease from its primordial value during “violent 
relaxation” in the process of galaxy formation and evolution}^, [95l-98|. The coarse grained 
phase space density for the species Uh is dehned asj”" 


Vh = 


Wt) 




3/2 


(III.24) 


The discrepancy with the result in ref. 4^ is due to the difference between = 0.12 and = 


0.14 


47|. 
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where the average is with the distribution function fh{qc) and Mhif) is given by flIII.3j) . 
Therefore the primordial coarse grained phase space density is given by 


Vh = 


Ur ^yy^fhhj)] 


5/2 


When the DM particles become non relativistic it becomes 

P^h _ ^ P^h 


Vh = 


MiiKl) 


3/2 


33/2 Af* < 


For a primordial phase space density, Vh, imposing the bound Vtoday < TPh 
constraint 


(III.25) 


(III.26) 
gives us the 


- 33/2M^ 
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^ DM 
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(III.27) 


If there is only one species the values of and cXj,^ can be inferred from the kinematics 
of dwarf spheroidal galaxiesand this constraint can be used to obtain a bound which 
complements those from CMB observations. 


A. Non-LTE freeze-out from multiple production channels: 


The distribution function fh{y) is a result of all the possible production channels that 
are kinematically allowed as the total gain and loss terms determine the solution of the 
cosmological quantum kinetic equation. Furthermore, each production process of a heavy 
neutrino species Uh may actually be the result of the decay of a heavy standard model 
particle (such as W, Z,t, fi, ■ ■ ■) into different channels and each channel may yield a different 
distribution function because of the kinematics, we refer to this as “kinematic entanglement”. 
Thus even for a single species Uh of heavy neutrino, its frozen distribution function may be a 
mixture of several contributions some colder than others as a consequence of the kinematics. 

While the general solution flll.36p is non-linear in the gain and loss terms of each channel 
because of the exponential terms in the relaxation rate 7 , if the condition flll.dOp is fulhlled, 
these exponentials can be neglected and the result is given by flll.dip which is linear in 
the gain rates allowing for an identihcation of the distribution functions associated with 
each production channel. In each production channel the heavy neutrino is kinematically 
entangled with a lepton and the non-LTE distribution function at freeze out will reveal this 
“entanglement” for example in the form of production thresholds. This aspect is another 
manifestation of “mixed DM”, in the sense that even for a given heavy species Uh, its 
distribution function at freeze-out is a result of different contributions from different channels 
each with different kinematics. The concentration of each component depends, among other 
factors, of the ratios of mixing angles for the different channels. If freeze-out occurs out 
of LTE, the distribution functions will maintain memory of the processes that led to the 
production and the kinematic entanglement, in LTE this memory is erased as the distribution 
function becomes the Fermi-Dirac distribution regardless of the production process. 

Separating the contribution from the different channels will also allow an assessment of 
the “coldness” of the heavy neutrino as a result of the particular production channel. An 
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explicit example of these phenomena will be studied in detail within the context of heavy 
neutrinos produced from pion decay in the next section. 


B. Summary of cosmological constraints: 


We now summarize the main cosmological constraints in terms of the decoupled distri¬ 
bution function, degrees of freedom and mass of a particular heavy species Once the 
distribution function has been obtained from the solution of the quantum kinetic equations 
these constraints inform the feasibility of such species as a suitable DM candidate. In the 
discussion below fh{y) is the total distribution function solution of flll.36p after freeze-out. 


Abundance: the fraction of DM from a particular species Uh must obey, < 1 
leading to the abundance constraint 






7.4 eD \gdh'' Jo 


dyy^fhiy) < 1, 


(III.28) 


Phase space density (Tremaine-Gunn): the result that the phase _space density 
only diminishes during “violent relaxation” in gravitational collapse j^. losl-^ leading 
to (IIII.27D yields the constraint 
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The most DM dominated dwarf spheroidal galaxies feature phase space densities within 
a wide range (see 98| and references therein) 
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Stability: To be a DM candidate the candidate particle must feature a lifetime larger 
than the Hubble time 1/Ho, namely tHq > 1. Heavy neutrinos feature various leptonic 
and radiative decay channels analyzed in section (JITl)- A conservative bound on the 
lifetime from the dominant leptonic decay is (see eqns. flH.521IH.59|) l 
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in particular for heavy neutrinos with Mh < 1 MeV the decay channel with the largest 
branching ratio corresponds to the “invisible” decay into three light active-like neutri¬ 
nos. This decay mode provides interesting possibilities even when the heavier neutrino 
is unstable, this will be discussed in section (jV]) below. 

Caveats: If indeed DM is composed of several species including a hierarchy of different 
masses, the individual free streaming lengths and phase space densities of a species Uh of 
mass Mh may not be the proper characterization. Although the different species do not 
interact directly with each other (or do so extremely feebly), they interact indirectly via 
the common gravitational potential which is sourced by all species. While a previous study 
addressed the gravitational clustering properties of mixed dark matter that study did not 
include cosmological expansion and should be deemed, at best, as of preliminary guidance. 
As far as we know, there has not yet been a consistent study of free streaming and phase 
space dynamics for mixed DM including cosmological expansion during the different stages. 
In particular on the important question of what is the correct cut-off scale in the linear 
power spectrum in the case of various components. This entails the study of the linearized 
collisionless Boltzmann equation including cosmological expansion. While a numerical study 
implementing the public Boltzmann codes may yield insight, to the best of our knowledge 
an analytical study for arbitrary concentrations of the different components clarifying the 
contributions of each species to the effective cutoff scale is still lacking. Until such study 
emerges, we will consider the free-streaming and phase space density constraints of the 
individual species as indicative. 


IV. HEAVY NEUTRINO PRODUCTION FROM PION DECAY. 

The previous section discusses in general terms the quantum kinetic approach to produc¬ 
tion and freeze-out of heavy neutrino species and the various processes that may produce 
them in the early Universe. It is clear from the discussion in section dlTl) that standard 
model interaction vertices that lead to the production of active neutrinos will also lead to 
the production of the heavy neutrinos as long as the processes are kinematically allowed. 
This entails a far broader range of production mechanisms than those that had been the 
focus in the literature and suggests that a hrm assessment of a particular candidate, such as 
a keV “sterile” neutrino requires a deeper understanding of all the standard model processes 
that may lead to their production during the various cosmological stages. 

We recognized the necessity to include hnite temperature corrections to masses and in¬ 
teraction vertices to obtain a reliable description. In this section we implement this program 
with a dehnite example: the production of heavy neutrinos from pion decay shortly after 
the QCD ^ase transition (crossover) into the conhned phase. 

In ref. 12^ the authors proposed to study hadronic contributions to the production of 
sterile neutrinos by considering the self-energy corrections to the active (massless) neutrinos 
in terms of correlators of vector and axial-vector currents. While this is correct in principle, it 
is an impractical program: the conhned phase of QCD is strongly coupled and the description 
in terms of nearly free quarks is at best an uncontrolled simplihcation, casting doubts on 
the reliability of the conclusions in this reference. 

Instead, in this section we study the production from pion decay by relying on the well 
understood effective held theory description of weak interactions of pions, enhanced by 
the results of a systematic program that studied hnite temperature corrections to the pion 
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decay constant and mass implementing non-perturbative techniques fro m chiral p erturbation 


108 1. There are 


theory, linear and non-linear sigma models and lattice gauge theory |l00 - 
at least three important reasons that motivate this study: 1 ) it is a clear and relevant 
example of the quantum kinetic equation approach to production and freeze out that includes 
consistently hnite temperature corrections. In fact this case is similar to the production via 
the decay of W, Z vector bosons the main difference being the three polarizations of the 
latter and kinematic factors. 2) Pion decay surely contributed to the production of heavy 
neutrinos in the early Universe if such species of neutrinos do exist: the QCD transition to its 
conhned phase undoubtedly happened, and the consequent hadronization resulted in baryons 
and low lying mesons, pions being the lightest, are the most populated meson degrees of 
freedom near the QCD scale. 3) The two leptonic decay channels vr —)■ tt —)■ ei^h feature 

different kinematics and thresholds, in particular the (V-A) vector-axial coupling results 
in chiral suppression of the e channel for vanishingly small M^. Since the heavy neutrino 
produced in the decay is kinematically entangled with the emitted leptonjH^, we expect 
that the distribution functions from each channel will display differences as a manifestation 
of this kinematic entanglement, thus providing an explicit example of memory effects as a 
consequence of “kinematic entanglement” in the non-equilibrium distribution function and 
dependence on the particular production channel. 

At temperatures larger than the QCD Phase transition Tqcd ~ 155MeU lOSj, quarks and 
gluons are asymptotically free. Below this temperature, QCD bound states form on strong 
inter action time scales, the lightest being the pion. Recent lattice gauge theory calculations 
|l08 | suggest that the conhnement-deconhnement transition is not a sharp phase transition 
but a smooth yet rapid crossover at a temperature ~ 155MeV within a temperature range 
AT ± lOMeU. This occurs in the radiation dominated epoch at t ~ 10 fisecs within a time 
range At ~ 2 — 3 usees. This is much larger than the typical strong and electromagnetic 
interaction time scales ~ 10“^^ secs implying that pions that form shortly after the conhning 
cross-over are brought to LTE via strong, electromagnetic and weak interactions on time 
scales much shorter than At. After pions are produced they reach LTE via vr — vr scattering 
on strong interaction time scales, they decay into leptons and maintain LTE via detailed 
balance with the inverse process since charged leptons and active-like neutrinos are also in 
LTE. In conclusion, for T < 155MeU, pions are present in the plasma in thermal/chemical 
equilibrium due to pions interacting on strong interaction time scales while the 

crossover transition occurs on the order of 2 — 3 fisecs. 

If the pions (slowly) decay into heavy neutrinos Uh with very small mixing angles, detailed 
balance for vr ^ luh, h = 4,5 ■ ■ ■ will not be maintained if the heavy neutrinos are not in 
LTE. As the pion is the lowest lying bound state of QCD, it is a reasonable assumption that 
during the QCD conhnement-deconhnemet crossover pions will be the most dominantly 
produced mesonic bound state and, during this time, pions will remain in LTE with the 
light active-like neutrinos and charged leptons by detailed balance tt ^ la^^a- We focus on 
heavy neutrino Vh production from vr —)• luh which is suppressed by -C 1 with respect 

to the active neutrinos and does not maintain detailed balance. 

A full study of sterile neutrino production through vr decay in the early universe requires 
various hnite temperature corrections and a preliminary study which focused on the pro¬ 
duction of neutrinos in the keV mass range |44| has implemented the hrst step. This study 
yielded a suitable warm dark matter candidate with free streaming lengths on the order of 
several kpc whereas it is expected that heavier neutrinos will yield a colder dark matter 
species. 
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The issue of a heavy neutrino production {Mh > MeV) through tt decay has not been 
addressed and is the main focus of this section. Through the two possible channels vr —)• 
pz/^; TT —)■ ez/g pion decay offers a wide kinematic window to the production of heavy neutrinos 
and provides a natural mechanism to produce a mixed dark matter scenario provided that a 
hier archy of heavy neutrino species exist and their production is kinematically allowed 

Furthermore, “kinematic entanglement” will yield different distribution functions for 
the different channels, thus providing an example of a mixed distribution for the same 
species. 

Whereas pion decay vr —)■ fiu, is one of the most ubiquitous mechanisms to produce neu¬ 
trinos in many terrestrial experiments, including short baseline experiments, this production 
mechanism has not been fully addressed within the cosmological setting. While we do not 
claim that the processes studied below are more or less important than the processes de¬ 
scribed in the previous section or discussed elsewhere in the literature, this study leads to 
a clear example of the concepts and methods described in the previous section, including 
the hnite temperature corrections. Furthermore, this production mechanism may yield a 
substantial (if not the dominant) contribution to DM below the QCD scale. 

In studying the production of heavy neutrinos, we make the following assumptions for 
the quantum kinetic equation that governs the Uh population build up: 


A substantial b ody of work has inv estigated hnite temperature corrections to the pion 


decay constant 102 - 104 106 . 107l |. To leading order in chiral perturbation theory, the 


correction is given by 


fn ^ m = fm 1 - 


T{tY 


U{0) = 93MeV. 


(IV. 1) 


This result will be needed for the quantum kinetics equation as we consider production 
beginning at Tqqd ~ l55MeV. The se re sults are well established but this sub held 
remains an active area of investigation |ll0 | and, for precision studies, these results will 
need rehnement. The assumption is that there are no pions present in the plasma above 
Tqcd 155MeV and that hadronization happens instantaneously in comparison to 
the neutrino production time scales. 

Finite temperature corrections to the masses of both the pions and charged leptons 
occur in th e pla sma. The pion mass including electromagnetic corrections is shown in 
hgure 2 of [lO^ . It is shown that between the temperatures of 10 and 150 MeV, the 
pion mass varies between 140 and 144 MeV. This change is relatively small and will 
be neglected for the calculation and we will use an average = 142MeV. Finite 
temperature corrections to the charged lepton mass are of (!VeT)[84-84. 86] and for 
muons this correction is only a fraction of its mass so it can be neglected. However, 
mass corrections to the electron may be substantial, but in this case we are interested 
in a kinematic window of large mass for i>h- For the purposes of this work, the charged 
lepton mass will be taken as a constant for both muons and electrons. The effects of 
charged lepton mass corrections will be investigated elsewhere. 


It is argued in [ill ] that the chemical potentials (including pions) are on the order 
of ~ 10~^eV for the temperature range we consider here. We assume that the lepton 
asymmetry is small and will be neglected in our calculations consistently with the 
neglect of chemical potentials in the quantum kinetic equations of the previous section. 























The mixing angle between active-sterile neutrinos in the presence of a matter potential 
will develop a temperature and lepton asymmetry dependence as shown in refs. [12, 
80, For vanishing lepton asymmetri es ( chemical potentials) and T/Mw -C 1 
there are no in medium MSW resonances [^. 1^ . In absence of lepton asymmetry 

the finite temperature in-medium potential lead to the following modification of the 
mixing angles (we consider mixing between one active-like and a heavy neutrino only) 


sin(26'^) = 


sin 29 


[sin^ 29 + (cos 29 + / A)^] 

where A = 5M^/2E = {Ml - Ml)/2E, and ~ ISO 

concerned with the temperature range E T 100 Meld, we hnd 

T f keV^ 2 
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We will be 


wva ~ 10-2 (—-— 

' VlOOMeV/ \mJ 
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The focus of this study is mainly on neutrinos with > few keV and sin(20) -C 1 
at temperatures below T ~ ISOMeld so the temperature dependence of the mixing 
can be ignored. In any event, the temperature correction to the mixing angles may, 
at most, lead to a slight quantitative change but not to a major modihcation of our 
main results, in particular the contributions from different channels with kinematic 
entanglement is a robust feature, as it will become clear from the analysis below. 

The effective low energy held theory that describes pion decay appended by hnite tem¬ 
perature corrections to the mass and decay constant yields the following interaction 
Hamiltonian (in the interaction picture). 


Hi 


T ^GFV„df,(T) J + h.c.] 
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where = i5o-7r(x, t) is the pseudoscalar pion current, Gp is the Fermi constant, Vud is the 
CKM matrix element, / 7 r(T) is the pion decay constant with hnite temperature corrections 
given by fllV.ip and the “havor” neutrino helds are related to the helds that create- 
annihilate neutrino mass eigenstates by the relation flll.ip . 

Pions decay into heavy neutrinos through the channels ; tt’*' —)■ e^Vh (and the 

charge conjugate). Because of the kinematics the heavy neutrino Vh is “entangled” with the 
charged lepton I in the sense that the gain and loss rates depend on the particular channel. 
Therefore for each channel we label the gain and loss rates and the distribution function 
with the labels Z, h. As discussed in detail in section ([ni) the total rates are the sum over all 
channels and the quantum kinetic equation inputs the total rates as per eqn. flll.34p 

The quantum kinetic equation for production is given by fill.331) with the total gain 
and loss rates summed over the kinematically allowed channels labeled by the corresponding 
charged lepton 1. The results for these rates are given in detail in appendix flAp . 
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and 


ril(?) = r 5 ( 9 ), (IV.6) 

where momenta and energies are replaced by their local expressions in the expanding cos¬ 
mology flII.30D and 


P±{t) = 


- Mf + Mif - 4MjM,;]V= ± M? + Mi 


2Ml 


2Ml 


and the constraint from the energy conserving delta fnnction 

K(p,t) = Ei(k,t) + E^(q,t) ; k=\p-if, 


(IV.7) 


(IV.8) 


As argued previously, the pious and charged leptons are assumed to be in LTE so 

that their distribntions will take the standard bosonic and fermionic forms 
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where Qc is the comoving momentnm and we neglect chemical potentials in this stndy. 

Inserting the distribntions into eqn. (1IV.5P and nsing Ei{k) = E.„{p) — Eh{q), the gain 
rate becomes 
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where are the solntions to the phase space constraints of eqn. (IA.17I1 and the time 
dependence of the physical momentnm has been snppressed. This integral is readily solved 
by a snbstitntion with the following result 
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(i—Eh{q,t)/T(t)^ETr/T(t) _j_ 


E-n=E^ iqf(t)) 


where 


Et(v(t)) = ^ 


£/.(<;/(*)) Mf + Mf - Mf) ± qs(t)X(M„ M,, Ms 


and the threshold function is dehned as 
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In the limit Mh —)■ 0 the bracket in the hrst line of fllV.llj) yields the usual factor — 

M^), which is the hallmark of pion decay vanishing in the Mi ^ 0 limit. 

In eqn. (lIV.lip . it is readily seen that r< depends on the ratio y = where Tq 

the temperature of the plasma today with the scale factor normalized today at {t = t^) with 
a(f*) = 1 and Qc is the comoving momentum. 

The threshold function is one of the signatures of “kinematic entanglement”; for hxed 
Ml, the threshold function vanishes at M^ = M^^ — Mi, for this value of Mh the two roots 
coalesce and the rate vanishes. This is important because for Mh ^ 36 MeV there are 
two production channels of heavy neutrinos, tt —)■ yi'h'M whereas for Mh > 36 MeV 

only TT evh is kinematically available. 

Restricting our study to the production process tt ^ Ivh and its inverse, we anticipate 
that the condition for production and freeze-out out of LTE (III.4011 will be fulhlled, this will 
be proven self-consistently below. In this case the solution of the quantum kinetic equation 
is given by flll.4111 . namely the quantum kinetic equation flll.33p simplihes to 

dnh{qc'-it) f ff\T ■[ A\ 

--= Motive] t), (IV. 14) 

where 

r£,(?c;i) = (iv.i5) 

l=fj,,e 

and the sum is over the kinematically allowed channels. In this approximation the linearity 
of flIV.1411 allows us to introduce a distribution function for each channel nih{qc] t) that obeys 

dnih{Qci T< / j.\ 

--= ^ihiQc] t), (IV.16) 

which will prove useful as it allows the opportunity to understand how the “kinematic 
entanglement” associated with each channel is manifest in the frozen distribution function. 
In this approximation the total distribution function is 

Uhiqc] t) nihiqc] t). (IV.17) 

l=fi,e 

and after freeze-out 

fh{qc) = Uhiqc] t = oo) = ^ fihiqc] t = oo), (IV.18) 

l=fj,,e 

in other words the total distribution function after freeze-out is a mixture of the contributions 
from the different channels. 

We emphasize that whereas the above dehnition is useful to learn the effects of the 
“kinematic entanglement” in the frozen distribution function, the cosmologically relevant 
quantities, such as coarse-grained phase space densities, free streaming length etc, are non¬ 
linear functions of the distribution function as they involve moments, therefore for these 
quantities only the total distribution function Uh fllV.lTp is relevant. 

It proves convenient to introduce a dimensionless parameter r that will play the role of 
time along with with a change of variables that introduces manifestly a factor of the Hubble 
factor (during radiation domination): 


?/(^) _ Qc ^ 

T(t) ~ To • ^ “ T(t) 


- = rH(t) ; H(t) 
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furthermore take the overall scale to be and dehne the following dimensionless ratios 


Mh Ml 

rrih = —— : mi = —— ; i = a.e . 


(IV.20) 


As argued previously, production from pion decay begins after the QCD crossover during 
the radiation dominated epoch. We will argue that freeze out occurs at temperatures T > 
lOMefd so that the Hubble factor is given by eg lIV.lQl for the entire period of production. 
With this change of variables and factoring the constants out of Eg. fllV.llj) leads to the 
dehnition of the overall scale of the problem: 
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where the dimensionless ratios mi^h are dehned in (IIV.20D . After the temperature falls below 
the hadronization temperature (T ~ 155MeH), the relativistic degrees of freedom {g{t)) 
remain constant until the temperature cools to below the /i mass (T ~ 106MeH) and again 
remains constant until cooling to below the electron mass. For 155MeV > T > 106MeH 
the relativistic degrees of freedom are q(t) ~ 14.25 and for 106MeV" ^ T > O.bMeV the 


degrees of freedom are g{t) ~ 10.75 [^. As will be argued, Uh production is complete by 
lOMeH and since g{t) has a small variation in the temperature range of production, g{t) is 
replaced by its average value g = 12.5. With this, the time variation of the overall scale Ah 
can be neglected, and using \Vud\ = 0.974[^ we hnd 


Aih ^ 0.13 


Hih\^ 

10-5 


2 I 2 

mi +mh 



(IV.22) 


The function 


C[mi, mh] = mf + ml- {m^ - mlf 


(IV.23) 


reveals the usual chiral suppression in the case of negligible neutrino masses as a conse- 
guence of the (V-A) coupling, being much larger for the /i-channel than the e-channel with 
C[me,mh]/C[mp,mh] ^ for mh ^ 0. 

The compilation of bounds on mixing angles of heavy neutrinos in ref. 112l | in combination 


with bounds from X-ray astrophysical data sl, 3^ suggests that <C 10 ® for fewkeV < 


Mh < 140 MeV conseguently Aih 1 for Mh in this range. 

Production from pion decay begins shortly after the QCD crossover at Tqcd — 155 MeV 
when pions form, therefore Mt,/Tqcd = tq < t < oo where tq — 0.92. 

Upon substituting the dimensionless variables, Hubble factor and overall scale the kinetic 
eguation may be written in the following form: 


1 dnih 
Aih dr 




T 



Ml \ 
drym ) 


\/y‘^ + mfj'^ ^Eh(qu)/T{t) _j_ 


In 



Eh{q,t)/T(t) ^ET,{p,t)/T(t) 



P=P+{t) 


P=P- (0 

(IV.24) 
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Upon evaluating the pion energy at the solutions of Eq. flA.lTp . in terms of the dimensionless 
variables dehned above flIV.19IIIV.20j) we obtain^ 


En{p±{t),t) 

T{t) 


2ml 


A 


Ih 


2 I 2 2 

y 


1/2 


±y5, 


Ih 


where we introduced 


(IV.25) 


Aih = 1 - + ml 


Sih = 


1 + mf + ml 


2m, 


2ml - 2mlml 


nl/2 


^Ih 


Ami 


nl/2 


(IV.26) 

(IV.27) 


From this equation, the rate can be integrated numerically to study the distribution 
function after freezeout and obtain the cosmological quantities discussed in section flllll) . 
However, before we proceed to a numerical integration of flIV.24|) we recognize important 
features of this equation that anticipate the behavior of the solution: 


• The prefactor, proportional to is small initially but the argument of the logarithm 
attains its largest value as the initial temperature T(to) > however, as r increases, 
the temperature decreases and the logarithmic term decreases eventually exponentially 
beating the growth of the prefactor. This analysis indicates that the production rate 
peaks as a function of time r and falls off fast. We conhrm this behavior numerically 
in Eg. (|3]). This hgure shows the fast rise and eventual fall off of the production rate 
which becomes exponentially suppressed for r > 10. This entails a freeze out of the 
distribution function: while the rate is exponentially suppressed, the total integral will 
remain hnite, thereby fulhlling the freeze-out condition (III.39p . In ref. ji^ a similar 
behavior was noticed for lighter neutrinos and we conhrmed numerically that a similar 
pattern holds for the whole kinematic range of in each lepton channel. 




FIG. 3: Production rate of a heavy Vh from vr —)• (eqn. (|IV.24h with I = y,e for M/, = IMeV. 


^ We do not include a label I in the result for Et^ to avoid cluttering of notation. 
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The expression flIV.25p with (IIV.26IIIV.27j) which is a consequence of the kinematics 
has very important implications on the “kinematic entanglement” of the heavy neu¬ 
trino. For fixed y, ruh the ratio flIV.25p is smaller in the p channel than in the e 
channel because both Ai^ and Sih are smaller, this implies that the rate is larger in 
the fi channel. This feature is also displayed in £g. (|2]). Therefore, we expect that 
nfj,h{y,T)/Afj^h > 'neh{y,'r)/Aeh and the distributions at freeze out to display this dif¬ 
ference. Furthermore, the function C[mi,mh\ flIV.23p in Aih (IIV.22P is also larger for 
the /x-channel than in the e-channel for small rrih- This is a consequence of the chiral 
suppression as a consequence oiV — A and is displayed in £g. (jiP . 




FIG. 4: for both channels in the kinematically allowed region of M^. 


The loss rate F^ is obtained from the gain rate (1IV.24P from the detailed balance condition 
(lIV.Sp and we can now proceed to integrate the quantum kinetic equation. However, we 
hrst conhrm that the gain-loss processes solely from pion decay (and recombination) lead to 
freeze-out out of LTE, namely we first confirm that the condition (111.401) is fulhlled. This 
is shown explicitly in figs. (|5]). These figures confirm that the condition (Ill.dOp is fulhlled 
provided Aih -C 1. The result (1IV.22D indic ates that this is the case for \Hih\^ <C 10 


1-5 


I n iuic aies Liiat ims is Lire case loi lu 

32-3^ and a compilation of accelerator bounds 112 


Cosmological bounds from X-ray data 
suggest that ^ 10“^ for fewkeV < < 140MeV, therefore the condition Aih ^ 1 

is satisfied guaranteeing self-consistently that the condition (Ill.dOp is fulhlled. 

In this case the population build up is obtained by integration of the rate eqn. (IIV.24P 
as per the general result (III.4ip . In order to exhibit clearly the contribution from 7r-decay 
we assume vanishing initial population, in this approximation any initial population must 
be added. For each channel we obtain the distribution function from direct integration of 
the gain rate 

nih{y,'r)= [ ^^(?/,r')dr'. (IV.28) 


'TO 


dr 


The results are shown in hgs (iGlITIIHp . Each hgure shows the distribution at diherent values 
of time (r) and, in all cases, the distribution has been frozen out at r ~ Tfr — 10 at which 
T{Tfr) ^ 14MeV, this is because the rate is being suppressed through the pion thermal 
distribution as displayed in hgs. ([3]). 
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FIG. 5: {1/kih) 7 ;i(y, r')fir' for the two channels . 


There are several features to observe from these plots, the hrst of which is that, for low 
mass neutrinos (M^ ^ IMeV), the distributions observed in are recovered. For heavier 
neutrinos (Mh ^ SOMeld) a very different behavior from the light species is observed. The 
light species features a vanishing distribution for small momentum (y 0), peaks at a 
particular momentum and falls off at large momentum whereas the heavier species has non¬ 
vanishing support at zero momentum (y = 0) and monotonically decreases as momentum 
increases. An intermediate behavior can be observed in the electron channel for a heavy 
neutrino with = 20MeV; the crossover in behavior occurs for — T(Tfr) ~ 10 — 

14MeV for T/r — 10 is the time scale at which the distribution functions freeze out, the 
freeze out time r ~ 10 is nearly independently of the value of M^. 

Furthermore, comparison of figs. ([6]) with those of figs. (IH]) is revealing. These two 
sets display the distributions from the electron and muon channels within the kinematically 
allowed window available for both channels, 0 < < 36 MeV, and it is in this comparison 

that the relevance of the distribution functions nih for I = fi,e become manifest. Two 
important aspects stand out: a) for small y the distribution function from the jj, channel 
is systematically larger than that for the e channel, this is a consequence of the kinematic 
factors Aih,6ih in flIV.25p . both are much smaller in the y, channel than in the e channel, 
this fact makes the energies in flIV.25p smaller in the y channel, therefore yield a larger 
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FIG. 6: Distribution function neh{y,T) from vr —y with = 1,10, 20,30MeF and vanishing 

chemical potentials. The solid line is the asymptotic frozen distribution. 


contribution since these are less thermally suppressed. This conhrms the analysis presented 
above, b) The distribution function for the e-channel has larger support for y > 1 than that 
for the /i channel, which is larger for y <1 but us strongly suppressed for y > 1 as compared 
to that of the e-channel. Both these features lead to the conclusion that the distribution 
function from the y channel yields a colder component than that of the e channel in the 
sense that the velocity dispersion obtained from is smaller than that obtained from Ugh, 
this is a manifestation of the kinematic entanglement and its consequences will be analyzed 
in detail below. 

The total distribution function is thus a mixture of a colder and a warmer component. 
Dehning the distribution function at freeze-out for each channel as 


fih{y) = nih{y,oo) 


(IV.29) 


where by r = cxo we mean r > Tjr — 10, the total distribution function is given by 


fh^y) 


ffih{y) 


0(36 Me V 


Mh) -|- Aeh 


fehiy) 

Ae/i 


0(141 MeV 


M,). 


(IV.30) 


The 0 functions describe the thresholds in each channel, the brackets 


flh/^Ih 


are actually 
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FIG. 7: Distribution function Uehiu, t) from vr —>■ with = 50, 75,100, ISOMeF and vanishing 

chemical potentials. The solid line is the asymptotic frozen distribution. 


the result of the numerical integrations as per eqn. flIV.24p and are given by the solid lines 
(r > 10) in hgs. ([MHD- 

For a given range of masses the only unknowns are the mixing matrix elements Hi^. 
Because the total distribution function is a result of the combination of several production 
channels, each one with possibly a different mixing matrix element Hih-, it is convenient to 
factor out one of these and rewrite the total distribution function flIV.30p in terms of ratios, 
namely 




10 


fh{y) = 0.13 

Jhiy) = C[m^,mh] 


-5 


fh{y) 


(IV.31) 


f^ih iy) 


A 




(\H, 


eh 




C[me,mh] 


0(36 MeY -Mh 
feh{y) 


A 


eh 


Q{UlMeV - Mh) 


(IV.32) 


where the coefficients C[mi,mh\ are given by flIV.23D and the brackets [//A] are the results 
obtained numerically and displayed by the solid lines in hgs. ([MH])- In this manner the various 
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FIG. 8: Distribution function r) from vr —)• fiUh with M/j = 1,10,20,35MeF and vanishing 

chemical potentials. The solid line is the asymptotic frozen distribution. 


constraints can be phrased in terms of an overall mixing matrix element and ratios, snch as 
which along with the kinematic factors C[mi,mh] determine the concentration 
of the warmer species (from the electron channel) in the mixture of colder and warmer 
distribntions from both channels. 


A. A Tale of Two Distributions: 

Closer examination of hgs. [61I7II8] reveal a transition in the shapes of distribntions as a 
function of nentrino mass. In the case of relatively light nentrinos with Mh < 1 —lOMeld, the 
distribntion fnnction vanishes at small momentnm, reaches a maximnm, and falls off with a 
long tail at high momentnm. The sitnation for heavier nentrinos, > SOMeC, produces 
a distribution function which features a non-vanishing plateau at low momentum which 
steadily falls off with increasing momentum. This latter distribution obviously produces a 
colder dark matter candidate as the distribution function has considerably more support at 
small momentum. The transition between these two types of distributions is controlled by 
increasing neutrino mass and an intermediate distribution can be seen in the intermediate 
mass range, lOMeV < Mh < SOMeV, illustrated by the distribution function in hgure [6] 
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for a neutrino mass Mh = 20MeV in the electron channel. For this Mh = 20MeV case, the 
distribution resembles a superposition of the two distinct distributions seen for the lighter 
{Mh < lOMeV) and heavier species {Mh > 30MeV). 

Unfortunately, only limited analytic progress can be made towards an understanding of 
the full distribution function but this proves enough to shed some ligM on what governs 
the transition between the different distribution shapes. In reference [^, it is shown that 
under appropriate approximations, the distribution function for light mass neutrinos takes 
the form 


nih{To,y) 


Mh<lMeV 


A 


1 + {-l)^+^e^y 
1 + ey 


exp{-ky/Aih) 

k 


xJfc(ro,|/) (IV.33) 


where 


Jk{ro,y) = 2ro 



( kAihT^ \ r y 

V ^yj \kAih) 


2y 

kAih 




(IV.34) 

where r(z/, z) is the incomplete gamma function and the details of this calculation are repro¬ 
duced in appendix IHl As detailed in appendix [Bl the main approximation employed towards 
producing this semi-analytic expression is Mh/M^^ 1. This approximation is equivalent 
to the assumption that the neutrino is produced ultra-relativistically {Eh/T ~ y) and this 
proves to be an excellent approximation for lighter species {Mh ^ IMeV). 




FIG. 9: Comparison of the exact distribution with a distribution producing ultra-relativistic and 
non-relativistic heavy neutrinos. 


The opposite type of approximation, where the neutrino is produced non-relativistically 
{Eh/T M/T + -^y'^), can be made but leads to a very unwieldy expression which is given 
by Eqs lB.4IIB.lll The production of non-relativistic heavy neutrinos leads to a natural source 
of non-trivial behavior at small momentum which explains the low momentum plateau. 
These arguments lend themselves to the interpretation that neutrinos with Mh < lOMeU 
are solely produced ultra-relativistically while some fraction are produced non-relativistically 
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as the neutrino mass increases. This argument is illustrated in hg. |9] where the distribution 
functions of IB.9IIB.lll which are results of ultra-relativistic and non-relativistic production 
approximations, are plotted against the full distribution. In £g. |9l the low momentum 
region is dominated by the non-relativistic result whereas the large momentum region hts 
quite well to the ultra-relativistic result. Both the ultra-relativistic and non-relativistic 
results fail outside the appropriate momentum regions: the non-relativistic approximation 
fails for higher momentum whereas the ultra-relativistic approximations fails at decreasing 
momentum. 

This analysis explains the origin of the features of the distribution functions obtained 
numerically in hgures O [TJ [H namely that the ultra-relativistic approximation serves as a 
good estimate for the whole distribution as the neutrino mass becomes negligible compared to 
the pion mass while the non-relativistic approximation serves to understand the appearance 
of the plateau at low momentum in the distribution. This means that for light mass neutrinos 
(M/i -C Mt,) all of the neutrinos that are produced are done so ultra-relativistically while 
the production of a more massive species leads to some fraction of neutrinos produced non- 
relativistically as well. The fraction of neutrinos which are produced relativistically increases 
as sterile neutrino mass decreases while the fraction which are produced non-relativistically 
increases as the sterile neutrino mass increases. 


B. Non-thermality 


The equation of state parameter, w(T), is given by eq IIII. Ill where w = 1/3, 0 correspond 
to ultra-relativistic and non-relativistic species respectively. Depending on the temperature 
of production and decoupling, a heavy neutrino could be ultra-relativistic, non-relativistic 
(M -C T , M S> T) or somewhere in between as freeze-out occurs and w{T) determines 
when a particular species becomes non-relativistic. 

The correct equation of state fllll.llD and velocity dispersion flIII.12p must be obtained 
from the total distribution function flIV.30p . because these are moments of the distribution 
function they are not simply the addition of the two components. 

However, it proves illuminating to define an equation of state for each channel 


wi{T) 


V 

P 


Jdy j ""^2 

1 

^ Jdyy^\jy‘^ + ^ fih (Qc) 


(IV.35) 


these serve as proxies to quantify the “coldness” of the species produced by the particular 
channel: from the discussion in section fillip . ^^w{T) is a generalization of the “adiabatic 
speed of sound” for collisionless DM, and in the non-relativistic limit w{T) —>■ (l/^)/3. 

The distribution function of a thermalized heavy neutrino would be given by the standard 
Fermi-Dirac distribution: 


fLTE{y]T) — 


(IV.36) 


^ I ’ 

and we compare the equation of state wi{T) (1IV.35P with that obtained from (]IV.36p . this 
comparison quantihes the “non-thermality” of the distribution functions from each produc¬ 
tion channel. 
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FIG. 10: wi{T) for both channels tt —>■ ^ ei^h compared to thermal distribution for 

Mh = 1,10,20, SOMeV in the kinematic window in which both channels are available. 


The equation of state wi{T) for each channel and for the thermal distribution are com¬ 
pared in hgs (llOmip as a function of 1/T for a range of different masses. These figures 
clearly display the non-thermality of the heavy neutrino species, furthermore, as anticipated 
by the discussion above, the distribution function from vr —)■ /izz/j, namely f^h yields a colder 
component than that from the e channel, a direct consequence of the “kinematic entan¬ 
glement” leading to a larger amplitude at small y for and both components are colder 
and becoming non-relativistic w{T) -C 1/3 much sooner than the thermal case. 

The total equation of state is given by the full distribution function (IIV.31IIIV.32j) . namely 


V 

w{T) = - 
P 


Jdy i^\, fh{qc) 

1 _ 

^ Jdy y^ \/y^ + ^ fh{qc) 


(IV.37) 


which depends on the ratio as this ratio varies between 0 and oo, it follows 

that w{T) for 0 < Mh < 36MeV interpolates between the two dashed lines corresponding 
to the muon and electron channels in hgs. (nnii. this is shown in £g. flT^ . For the mass 
range 141 MeV > Mh > 36 MeV only the electron channel contributes and w{T) is given by 
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FIG. 11: wi{T) for vr —)• ei^h compared to thermal distribution for = 50, 75,100,130 MeV. 


the results displayed in fig. llTTll . 


C. Cosmological constraints: 


Having obtained the distribution function at freeze-out, we can now implement the general 
results obtained in section fllllj) and establish the allowed regions within which the various 
cosmological constraints discussed in section fillip are satished. 


Abundance: with fh{y) and fhiy) given by flIV.321 IIV.3ip respectively and taking 
= 2 the abundance constraint (1III.28P becomes 


2.81 



/ y^fh{y)dy<l. 


(IV.38) 


This is a function of the ratio \Hf.h\^/\H^h\^ and also of through the coefficients 
C[mi,myl\. This bound was calculated in ^ for < IMeV this mass region is 
completely dominated by the muon channel, whereas here we allow masses up to 
Mh < Mjr — Ml which translates to Mh < 142MeV for the electron channel and 
Mh < 36MeV for the muon channel. 
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FIG. 12: Equation of state with full distribution function for the ratios \Hf.h \^= 0.01,1,100 
as a function of {MeV)/T for = l,35MeV. w{T) interpolates between the results for /r and e 
channels. 


• Phase space (Tremaine-Gunn): 

Using the smallest phase space value of ref which comes from the Fornax dwarf 
spheroidal galaxy: {p/o'^)today = 2.56 x 10“^ (fceU)^, and using the general result 
(IIII.29P with = 2 we obtain the constraint 


Mh y 

keVy V 10-5 J 


> 0.0038 


I dy y^fhiy) 


3/2 


Jdyy^fhiy) 


5/2 ■ 


(IV.39) 


• Stability: The “conservative” stability constraint (1III.3211 is independent of the dis¬ 
tribution function. 


These bounds and allowed parameter space are shown in Eg. flT^ along with 
the parameters of heavy neutrinos from the recently reported X-ray signals jssl, [s^ 
(M/j = 7.1 keV, \Hah\‘^ = 7 x 10“^^). An important observation about this hgure is 
that the inclusion of both channels in the total distribution function leads to (marginal) 


consistency with the claimed X-ray signal. This differs from ref. which claimed 


that consistency with the X-ray data was supported in the muon channel but not the 
electron channel. The main differences between these results and those of ref ^ are 
that we generalize the results to arbitrary mass (as opposed to ;< IMeU) and we in¬ 
clude both production channels as opposed to considering each channel independently. 


An interesting aspect in the two top figures are “kinks” in the abundance and 
Tremaine-Gunn (phase space density) line, this is a consequence of the kinematic 
thresholds in the full distribution function flIV.30p . 


• Free streaming scale: the free streaming wavevector and length scale are given by 
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FIG. 13: The bounds on from abundance, stability and phase space constraints. The 

allowed regions determined from Eqs (liy.38 IV. 391 , IIII.32I) are shown and the parameters which 
potentially explain the 3.5 keV signal [33l. 341 are also shown. The kink is a consequence of the 
thresholds. The allowed parameter space is within the region bound by the three lines. 


eqns. (lIII.16IIIII.19p respectively, namely 


A^.,.,(0) ^ 5.3 



I dyy^fhjy) 

\ I dyy‘^fh{y) 


kpc, 


(IV.40) 


again this is a non-linear function of the mass and for a species Uh produced by a 
single channel it would be independent of the mixing angle, however if there are several 
channels, as is the case in tt decay, it depends on the ratio of mixing angles. Fig. ffT4|) 
displays as a function of for various ratios. 


For Mh < 1 MeV the p (coldest) channel dominates, and A/s(0) is insensitive to the 
ratio, however, as — 30 MeV the contribution of the e channel becomes substantial 
and dominates above the threshold at — 36 MeV, the kink in the hgure is a result 
of this threshold. Therefore A/s(0) interpolates between that of the y and electron 
channels as a function of the ratio, just as the equation of state interpolates between 
the colder (p) and warmer electron dominated components. 
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FIG. 14: \fs{Q)/{kpc) as a function of for the ratios = 0.01,1,100. Right panel 

zooms in to highlight the kink because of thresholds and interpolation. 


D. Other processes in the same temperature range 

In this section we focused on a detailed presentation of heavy neutrino production from 
pion decay after the QCD hadronization transition, primarily as a clear example where 
the hnite temperature corrections of the pion decay constant and mass have been previously 
studied in the literature. However, as the discussion in section ([11]) highlights, there are many 
processes that produce heavy neutrinos via charged and neutral current vertices provided 
the kinematics is favorable. In the temperature range just explored T ~ 150 MeV muons 
are thermally populated and muon decay fi is also a production mechanism that is 

available provided Mh is in the kinematic window for the three body decay. However, this 
process, is sub leading in the temperature range T < M^, this is because the ratio of decay 
rates 



however while pions are only available as a production channel below Tqcd, muons, on 
the other hand, are thermally populated at larger temperatures therefore they contribute 
substantially to the production of heavy neutrinos with masses within the kinematic window. 
Their contribution to the total abundance of heavy neutrinos merits a deeper study, along 
with all the other mechanisms described in section (P. Furthermore, light neutrinos are 
thermally populated in a much wider temperature range and t Uh is the three 

body “fusion” process described by the gain term flH.2ip that yields heavy neutrinos at 
temperatures T > M^. The inverse process, Uh —t is the decay process described 

by the loss term flH.22p . which contributes even at zero temperature and describes the decay 
of the heavy neutrino (III.55p . The detailed study of all of these processes clearly defines an 
extensive program to assess reliably the production of heavy neutrinos in cosmology. 
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E. Comparison with Dodelson-Widrow [15l| : 


Although the Dodelson-Widrow fPW) [l^ (non-resonant) mechanism of sterile neutrino 
production via active-sterile oscillations has been recently shown to be inconsistent with 
cosmological data at > 99% confidence level [i^, it is illustrative to compare the results 
obtained above for the non-equilibrium distribution function to the (DW) case for further 
understanding of its cosmological consequences. 

The (DW) distribution function is 


fowiy) — 


/S 


ey + 1' 


where f3 is determined by saturating the DM abundance, namely (l^. 


(3 — VtoMh^ 


94 eD 


ll.deD 


(IV.42) 


(IV.43) 


Mh ) Mh 

Furthermore, we will focus our comparison on the mass scale — 7 keV because this scale 
is of observational relevance jss], 341 and we consider the /i channel since for ~ 7 keV this 
channel features the largest branching ratio. The distribution function obtained from pion 
decay, fh{y) is defined by eqn. fHV.29jl . An important ingredient in the comparison are the 
following results: 


1 


1.803 

(IV.44) 

1.830 

(IV.45) 


7^^ih J 

therefore the total integral of the distribution functions divided by their prefactors is ap¬ 
proximately the same, however as gleaned from fh{y)/A^h is strongly peaked at small 

momenta, and as discussed in the text the distribution function is fairly insensitive to Mh 
for Mh < IMeV. We provide two different manners to compare the distributions; i) fixing 
[3 and A^h to saturate the DM density in both cases ii) for fixed values of Mh and mixing 
angles \H^h\‘^ extracted from the analysis of ref. 42| to obtain A^h but consistent with the 
(DW) scenario we keep the value of (3 that saturates D_dm- 

Comparison 1: both (3 and A^/^ are fixed to yield the total DM density, namely the 
fraction fllll.lOD is = 1 in both cases. This yields 


f3 = 1.611 X 10 


-3 


Af,h = 1.600 X 10 


-3 


(IV.46) 


with these values we display both distribution functions in figJT^ 

Although the integrals of the distribution functions are the same, obviously fh{y) is 
sharply localized at smaller momenta, therefore yielding a colder distribution. 

Comparison2: for this case we keep (3 = 1.61 x 10“^ so that Jdw saturates the abun¬ 
dance bound, but A^h is now extracted from the upper bounds on the confidence band for 
the (DW) species of fig. (4) in ref.ji^, identifying the mixing matrix element with 

sin^(26') in this reference. We read from this fig. the values 

Mh-lkeV ; ~ lO"" (IV.47) 
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FIG. 15: The distribution functions fowiv) (dashed line) and fhiv) solid line, for j3 K. In this 
case both distribution functions are fixed to saturate the DM density. = 7 keV 


with these values we find 

A^h = 3.18 X 10"^ (IV.48) 

the comparison between the distribution functions is displayed in hgJTbl In this case the 
distribution function fh{y) yields a fraction = 0.205 to the DM abundance. 



FIG. 16: The distribution functions fowiv) (dashed line) and fh{y) solid line, for j3 = 1.61 x 
10“^ ; = 3.18 X 10“^. In this case the (DW) distribution function is fixed to saturate the 

DM abundance, whereas A^/^; Mh = 7 keV are fixed by the upper limits of the conhdence band for 
(DW) in hg. (4) of ref. 4^. For this case fh yields a fraction T = 0.205 of the DM density. 


Again it is clear that even when the distribution function yields a smaller fraction, it 
features a larger contribution at smaller momenta and falls off faster at larger momenta, 
this feature makes this species colder than (DW) with an abundance that while smaller 
than (DW) is fairly substantial. 

A further illuminating comparison is obtained from the free streaming length fllll.lQj) . 
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which is independent of the normalization factors /3,Af^h respectively. We find 

^^w)(o) _ 2.7116 fcpc (IV.49) 

Aj^(O) = 1.0231 fcpc, (IV.50) 

the low momentum enhancement of fh{y) yields a much colder distribution with a much 
shorter free streaming length as compared to the (DW) case. 

We conclude that while the (DW) mechanism seems to be ruled out as a sole production 
channel of sterile neutrinos, the comparison with the non-equilibrium function obtained from 
pion decay in the dominant channel offers a useful insight into the properties of sterile neu¬ 
trinos produced via this mechanism and suggests that the abundance from this alternative 
scenario could be a substantial contribution to the DM component. 


F. Comments: 


We have implemented the results on the finite temperature dependence of the pion de¬ 
cay constant and pion mass from a substantial body of work on chiral perturbation theory 
including reso nances and linear a nd non-linear sigma models including contributions from 
vector mesons 102Ml04 IIOGL Il07l |. Taken together this body of results offer a consistent 


description of the temperature dependence. However, we recognize that there remain uncer¬ 
tainties inherent in an effective description of hadronic degrees of freedom b ut qu estioning 
the validity of ChPT is beyond the scope of this wor k. Recently a latti ce study |ll0j | reported 
results some of which are at odds with those of refs. 102l - 104 . 106l 107l | such as an increase in 
the pion decay constant and a decrease of pion temperature near the QCD crossover scale. 
While a confirmation or rebuttal of these results and/or a resolution of the controversy is 
awaiting, we can speculate on the impact of the results of this reference, if these hold up. 
First: if the pion decay constant is larger this obviously results in a larger production rate, 
secondly: if the pion mass is smaller, this also leads to a larger production rate (less thermal 
suppression of pions in the medium) and crucially to a delayed freeze out with a longer stage 
of production as pions remain populated in the medium for a longer time. Lastly: if the 
pion mass falls below the muon channel threshold at high temperature, the production is 
solely through the electron channel resulting in a warmer distribution, as the pion mass in¬ 
creases towards smaller temperature, the muon channel opens up and the muon production 
channel with a colder component becomes available thus confirming the argument on the 
“mixed” contributions to the distribution function. All of these aspects bolster the case for 
consideration of pions as an important production mechanism, and the last point in partic¬ 
ular, bolsters the argument on kinematic entanglement. We thus conclude that the features 
imprinted on the non-equilibrium distribution function from the various decay channels and 
mixed nature of the distribution function is a robust qualitative prediction. 


V. UNSTABLE HEAVY NEUTRINOS: CASCADE DECAYS INTO STABLE DM 

Unstable heavy neutrinos with lifetimes much smaller than 1/Ho do not play a direct role 
as a viable DM candidate, however they can decay via a cascade into lighter stable heavy 
neutrinos that could be viable candidates. To discuss this scenario in more concrete terms, 
let us consider a hierarchy of two heavy neutrinos Uhi, with Mh^ S> and assume that 
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the heavier, is in the kinematically allowed window that allows its production on-shell 
from TT decay, namely tt —)■ Ivhx- If so after being produced, Vhi will cascade decay into 
Vh 2 + leptons on a time scale ~ , namely the intermediate heavier Uhi yields yet another 

production channel to the lighter 

TT —)■ li'h2hh , (V.l) 


where li,l 2 are other charged or neutral leptons. In this process the intermediate goes 


on shell and the cascade is mediated by 
of the stable(r) species i'h 2 is given by [6 


resonant decay. In this scenario the production rate 

la 


^ T^hihh) , (V.2) 

where Br^Uhi —t i'h 2 hh) is the branching ratio. If decays into a lighter heavy 

it can also decay into the active-like- light neutrinos the decay amplitude for 

^mi^m 2 ^m 3 HU, wliereas the amplitude to decay into Uhi —t i'h 2 hh oc therefore Br oc 
and oc namely is suppressed by an extra factor |hrp. If the lifetime of Uhi is 

> 10^ secs it can decay into active-like neutrinos well after Big Bang Nucleosynthesis (BBN) 
avoiding the constraints on the number of relativistic active neutrinos during BBN providing 
a late injection of neutrinos into the cosmic neutrino background well after BBN. A lifetime 

> 10^^ secs would inject a lighter heavy neutrino as a DM candidate after matter radiation 

equality, just when density perturbations begin to grow under gravitational collapse. Two 
specihc examples illustrate these possibilities: a) consider Mh < 1 MeV from the discussion 
in section (III Ell the decay channel with largest branching ratio (~ 99%) is the “invisible” 
channel —)■ (see eqn. (111.5511 1 with a lifetime 


'j- 


10 ^ fMeV\\ 

\ W) 


(V.3) 


with Mh < 1 MeV ; ^ 10“® the heavy neutrino decays into active-like neutrinos 

after matter-radiation equality “injecting” a non-LTE component in the cosmic neutrino 
background after neutrino decoupling, b) The decay rate into a lighter is suppressed by 
another power of therefore a Uhi with Mh^ < 10 MeV ; \Hmh\‘^ ^ 10“® can decay into 

Z//J 2 with Mhi — few keV also after matter-radiation equality, now providing a heavy neutrino 
with a keV mass range as a DM candidate just at the time when density perturbations begin 
to grow. This latter case may yield an excess of positrons, however, to assess whether this 
is observationally signihcant requires a deeper study. 

Obviously these are conjectures that merit a far deeper analysis, however this scenario is 
similar to that posited in ref.[^ of a hierarchy of heavy degrees of freedom decaying in a 
cascade into lighter species that may act as stable(r) DM candidates. 


VI. SUMMARY, DISCUSSION AND FURTHER QUESTIONS 

The main premise of our study is that if sterile neutrinos are a suitable extension beyond 
the Standard Model in which these mix with active neutrinos via an off-diagonal mass matrix, 
diagonalization of the mass matrix to the mass basis implies that heavy neutrinos mass 
eigenstates couple to standard model leptons via charged and neutral current interactions. 
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The same processes that produce active-like light neutrinos also produce heavier neutrinos 
if kinematically allowed, albeit with a much smaller branching ratio determined by small 
mixing angles. We study the production of heavy neutrinos via standard model charged 
and neutral current interactions under a minimal set of assumptions: small mixing angles 
with flavor neutrinos, standard model particles are in local thermodynamic equilibrium. We 
obtain the quantum kinetic equations that describe their production to leading order in the 
small mixing angles and give the general solution in terms of gain and loss rates that obey 
detailed balance. A wide range of charged and neutral current processes available throughout 
the thermal history of the Universe lead to cosmological production of heavy neutrinos 
including the possibility of production from collective excitations and “rare” processes in 
the medium such as plasmon decay, and “inverse” processes such as 71 /^ 

We discuss the general conditions for thermalization and argue that heavy neutrinos 
with lifetimes > 1/Hq (the Hubble time scale) freeze-out with non-equilibrium distribution 
functions. We generalize the concept of mixed DM to the case in which a single species of 
heavy neutrinos is produced by different channels and argue that in each channel the heavy 
neutrinos produced are kinematically entangled with the lepton produced in the reaction. 
If the distribution function freezes ont of local thermal equilibrium it maintains memory 
of this kinematic entanglement in the form of a colder or warmer distribution fnnction as 
compared to other channels. We quantify the “coldness” by obtaining the equation of state 
parameter for each channel, which serves as a “proxy” for the velocity dispersion of the DM 
particle when it becomes non-relativistic. If several channels contribute to the production of 
a particular species, the total distribution function is a mixture with components produced 
from the different channels, these may be colder or warmer as a conseqnence of the kinematic 
entanglement. The concentration of each component depends on the kinematics and the ratio 
of mixing angles in each channel. 

We summarize the abundance, phase space density and stability constraints that a suit¬ 
able DM candidate must fulhll in terms of the total distribution function at freeze-out 
and discuss clustering properties such as primordial velocity dispersions and free streaming 
lengths. We compared the quantum kinetic framework with other treatments available in 
the literature, recognizing the necessity for a consistent non-perturbative treatment in the 
case of possible MSW resonances in the medium. 

An important conclusion of this analysis is that, whereas many efforts focus on some 
temperature scale and on particular production processes, in order to reliably assess the 
feasibility of a particular heavy neutrino DM candidate, all possible production channels 
of this species must be carefully analyzed throughont the thermal history of the Universe. 
An immediate consequence of this principle is to alter the “standard” production mecha¬ 
nisms - Dodelson-Widrow, Shi-Fuller, scalar decays - all of which assume a vanishing initial 
population. The example of pion decay provides a unambiguous source of sterile neutrino 
population which will modify the standard results in a prescription described in . Further 
work is in progress on this front. 

We argued that the final distribntion function of heavy neutrinos after freeze-ont is indeed 
a mixture of the various contributions to it from the different production processes that are 
kinematically available. Only in the case of a fully thermalized population will the “memory” 
from the different processes be erased, but non-equilibrium distribution functions will have 
imprinted in them the kinematic entanglement from the different processes. 

As an explicit example we studied the production of heavy neutrinos from charged pion 
decay after the QCD crossover into the hadronized phase within the effective held theory of 
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weak interactions of charged pions, including finite temperature corrections to the pion decay 
constant and mass. Pion decay is one of the main sources of neutrino beams in accelerator 
experiments and pions being the lightest hadrons formed after the QCD crossover, their 
decay surely contributed to the cosmological production of heavy neutrinos. While not 
claiming that this process is more or less important than others (a claim that requires 
a detailed assessment of the other production mechanisms) it provides a wide kinematic 
window M/j <140 MeV from two different channels and offers a clear example of kinematic 


entanglement: the distribution function from the /i channel is distinctly colder than that 
of the electron channel and the total distribution is a mixture of both. We obtain the 
allowed region of parameters that fulhll the abundance, phase space density and stability 
constraints, these are displayed in hgs. flT^ for various values of the ratio \Heh\‘^/\Hfih\‘^, 
the boundaries of these regions reveal the thresholds for the different channels, a hallmark 
of kinematic entanglement. The equation of state (or alternatively velocity dispersion) and 
free streaming length interpolate between the colder component and the warmer component 
from the fr and e channels respectively, as a function of Mh and the ratio , 

clearly displaying the mixed nature of the total distribution function. 

We conjecture that heavy neutrinos with lifetimes <C 1/Hq may decay into active-like 
neutrinos after neutrino decoupling injecting light neutrinos out of equilibrium into the 
cosmic neutrino background, and for Mh < 10MeV and \Hih\^ < 10“® may decay after 
matter radiation equality into another heavy but lighter and stabler neutrino that may be 
a suitable DM candidate. 

Furthermore, we have provided a detailed comparison between the non-equilibrium (non- 
thermal) distribution function obtained from pion decay and that of the Dodelson-Widrow 
scenario. The comparison was carried out within two different scenarios: i) saturation of 
DM abundance for each case, and ii) the parameters Mh, \Hih\ were extracted from the 
upper bounds of the conhdence band in the analysis of ref. ^ . In the second case the 
distribution function from pion decay yields a substantial fraction of DM, in both cases 
the distribution function from pion decay yields much shorter free streaming lengths. This 
comparison suggests that production via pion decay could yield a substantial contribution to 
the DM abundance with a species that is colder than a thermal species but without invoking 
resonant production via a leptonic asymmetry. 

The wide range of production mechanisms available throughout the thermal history of 
the Universe explored in this study suggests the necessity of an exhaustive program to assess 
their contributions to DM from heavy neutrinos. 


Further Questions. 

Our study raises important questions that merit further and deeper investigation: i:) we 
recognize that the possibility of MSW resonances in the medium would require going beyond 
the leading order in the mixing matrix element to obtain the quantum kinetic equations (this 
is also a caveat of the approach in ref.|25j), one possible avenue would be to obtain the non¬ 
equilibrium effective action for the neutrino sector by tracing over the remaining degrees of 
freedom of the standard model (quarks, charged leptons, vector bosons) assumed to be in 
LTE. ii:) we argued that collective excitations in the medium could lead to the production 
of heavy neutrinos at high temperature, for example via plasmon decay, this is a cooling 
mechanism of stars in advanced stages of stellar evolution, and its high temperature coun¬ 
terpart in the early Universe may be a suitable production mechanism of heavy neutrinos. 
This is an intriguing pos sibility that requires a thorough analysis implementing the hard 


thermal loop program [83|, l8J, l8g to obtain the plasmon dispersion relations and couplings 
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to neutrinos, iii:) if there is a hierarchy of heavy neutrinos there is the possibility of mixed 
DM: the contribution from different heavy neutrinos providing cold and warm components 
depending on their masses and distribution functions. In this case (unlike the case where 
mixed DM arises solely from one component but from various production channels) it is far 
from clear what is the effective free streaming length and coarse grained phase space density. 
Since the free streaming length determines the cutoff scale in the linear power spectrum of 
density perturbations it is important to obtain a reliable assessment of its interpretation in 
the case of mixed DM, is it possible that a very heavy neutrino (cold DM) and a lighter one 
(hot DM) can mimic warm DM?, is it possible that such mixture would lead to the same 
power spectrum as one single species of mass ~ 7keV?. In order to shed light on these 
questions, the collisionless Boltzmann equation for several DM components in an expanding 
cosmology must be studied. Similar questions apply to the effective coarse grained phase 
space density, as this quantity is observationally accessible (or inferred) from the kinematics 
of dwarf spheroidal galaxies and provides a powerful constraint on a DM candidate inde¬ 
pendently of abundance. We expect to report on some answers to these questions in later 
studies. 
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Appendix A: Quantum Kinetic Equation 

In this appendix we set up the quantum kinetic equation describing the sterile neutrino 
pop ulation arising from the reaction tt^ ^ l^i^si^s)- The kinetics were originally set up in 
|44| but is repeated here for completeness. This process will occur when the plasma is at 
temperatures below the QGD phase transition and it is acceptable to use an effective field 
theory which treat pions as the fundamental degrees of freedom. The low-energy effective 
interaction Hamiltonian responsible for this process is given by eqn. flIV.4D including hnite 
temperature corrections to /^. 

The population buildup of the heavy neutrinos is described with a quantum kinetic equa¬ 
tion given by form fill. 81) where the gain and loss terms are calculated by writing the quantum 
mechanical transition amplitudes from an initial state z to a hnal state /, and in¬ 

tegrating over kinematic region of phase space. With the effective Hamiltonian flIV.4D . the 
relevant reactions are displayed in Eg [T71 

The gain terms are due to the decay process li'i which starts from an initial state 

with n^(|7) pion quanta and na{k), nh{q) charged lepton and heavy neutrino mass eigenstate 
respectively while the final state has — l,na{k)+ l,nh{q)+ l quanta for each respective 
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FIG. 17: The gain/loss terms for the quantum kinetic equation describing vr'*' —)• 


species. The Fock states for the decay process are given by 

\i) = \nnip),ni{k),nh{q)) ; |/) = \n^{p) - 1, nj{k) + 1, nh{q) + 1) . (A.l) 

In a similar fashion, the loss terms are obtained from the recombination lui —)• which has 
an initial state with n^(p), ni(fc), n/i(g) pion, charged lepton and nentrino qnanta respectively. 
The hnal state is popnlated with nT,{p) + l^ni{k) — l,nft(g) — 1 of the appropriate qnanta 
and the Fock states for the loss process are given by 

N) = \n^{p),niik),nh{q)) ; |/) = \n^{p) + 1, ni{k) - 1, nh{q) - 1). (A.2) 

Where the relation between the nentrino mass and flavor eigenstates are given by fill. ID : 

A standard calcnlation yields the transition amplitndes at tree-level and the matrix ele¬ 
ments for the gain term (decay) is given by ( "H/ is the interaction Hamiltonian density) 




fi I gain 


—i 


d^x{n.„{p) - l,ni{k) -b l,nh{q) + l\Hi{x)\n^{p),ni{k),nh{q)) (A.3) 




271 U^f^{q,ai)fVLV\k,a2) 
,/8E^{p)Ei{k)E,{q) 


X - EM) - Ei{k)\/nM)^^ - ni{k)^/l - n^{q) 

where W, V are Dirac spinors and ai^2 are helicity qnantnm nnmbers, while the transition 
amplitude for the loss term (recombination) is given by 


Afilioss = -M d^xM+ip)+ l,nj{k)-l,nh{q)-l\'Hi{x)\nM),M^^)^MQ)) (A.4) 




Ih 


277 V\k,a2)fVLU^^{q,ai) 

7^ MEAp)Ei{k)E,{q) 


X - Ey{q) - Ei{k))^/nM) + l\/ni{k)^nh{q) ■ 


Summing over the hnal states leads to the transition probability per unit time {V is the 
quantization volume) 
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/c,p,cri,(J2 


gam 


\Afi 


^ (27r) 

‘2Eh{q) 

X n^(p)(l 


d^p 2\H^,^\^\Vud\^Gy^ 

(27r)3 2E^ip)2Eiik) 


Tr[fPL{<l + Mh)fVL{$ - Ml) 


ni{k)){l - n^(g))(5(E^(p) - Eh{q) - Ei{k)) (A.5) 


(27r) f d^p 2\H^i,\^\Vud\^Glf^ 


2En{q)J (27r)3 2E^{p)2Ei{k) 


2[2{p ■ q){p ■ k) -p'^iq-k)] 


X n^(p)(l - nz(/c))(l - n,,(g))(5(E^(p) - E^{q) - Ei{k)) 


(A.6) 


where T is the total interaction time, not to be confused with temperature, and k = \p — q\. 
The loss term is calculated in the same way with the replacement n-j^ —)■ l + n-,^ and \—n^n 
for leptons. With these steps and explicitly evaluating the energy/momentum conservation 
leads to the quantum kinetic equation governing the sterile neutrino population: 


where 


dnh{q-,t) _ dnh{q;t ) , dnh{q;t ) , 
dt ~ dt dt 


(A.7) 


dnh{q;t) 


2n 


d?p \Mfi 


dt 2Ehiq) J (27r)3 2E^{p)2Ei{k) 


nn{p){l - ni{k)){l - nh{q; t)) 


X S{E^{p) - Ei{k) - Eh{q)) ; k=\p-q], 


(A.8) 


dnh{q-,t) I 

dt 


27r f d^p |Ad/iP 

^ 2Ei,{q) J 2E^{p)2Ei{k) 

X S{E^{p) - Ei{k) - Ei,{q)) ■ k 


{1 + n^{p))ni{k)nh{q-,t) 

= |P - ^ , 


(A.9) 


and the averaged \Aifi\‘^ matrix element is given by 

= 4\H^,f\VudfGlf^ [2{p ■q){p-k)- p\q ■ k)] . 


Therefore 

= r<(g)(l -nfe(g;t)) - V>{q)nh{q]t ), 

with the gain and loss rates given by 


rX?) 


rXX 


27r f d^p |Ad/ip 

^ 2E,(g) J ^ 2E^{p)2Ei{k) 

X 5{E^{p) - Ei{k) - Eh{q)) 

271 r d?p |Ad/ip 

^ 2Eu{q) j ^ 2E^{p)2Ei{k) 
X S{E^{p) - Ei{k) - Eh{q)) ; k 


n^{p){l - ni{k)) 


(1 + n^{p))ni{k) 

= \p-^. 


(A.IO) 

(A.ll) 


(A.12) 


(A.13) 


Performing the angular integration using the delta function constraint we hud 
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T<{q) = 


\H^h\^\Vud\^Glf^ Ml{Mf + Ml) - {Ml - Ml)^ 


X 


Stt 

dpp 


Ip- 


qEhiq) 

n^{p){l-ni{\p-^)) 


(A. 14) 


r>(g) = 


2\2 


H^nl^lVudl^Cyi Ml{Ml + Ml) - {Mf - Ml) 


X 


Stt 

dpp 


qEh{q) 

{I + n^{p))n-i{\p - 


Jp_ ^/pETMl 

The integration limits p± are obtained from the constraint 

I(Ip 1 - k1)" + Mfp-‘ < EM - EM < [(IpI + k1)" + Mfp\ 

The solutions are given by 


(A.15) 


(A.16) 


P± = 


Eh{q) 


2Ml 


ml 


Ml + M^ 


2\2 






Mf + Mf] 


2Mf 


(A.17) 


Note that these bounds coalesce at threshold when Mf — Mf + Mf — 2MfMf and the 
population change vanishes as expected. 

Using the relations 

1 + ?^7r(p) = ; 1 — ni{k) = ni{k) (A.18) 


and using the energy delta function constraints, the detailed balance condition follows, 
namely 

r<(g)e'®U9)/T _ r>(g). (A.19) 

These results are extended to cosmology by replacing the momentum with the physical 
momentum, q ^ Qf = qc/ci{t) where physical energy and momentum is measured by an 
observer at rest with respect to the expanding spacetime. 


Appendix B: Approximate Distributions. 

The exact rate equation IIV.24I is unwieldy and only limited analytical progress can be 
made. Recasting the rate equation in terms of the dehnitions of IIV.20IIIV.26IIIV.27I and a 
dimensionless variable for the neutrino energy 


£ = Eh/T 


A/f2 ^_ 

+ 


hy2 


(B.l) 


leads to the following form of the rate equation which is more suitable to manipulations: 
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■i^{5ihy + /Aihs) 
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1 + exp (-i^{5ihy + {Aih - 2ml)e) 
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2^ 


-dihy + Aihs) 


1 + exp 




-6ihy + {Aih - 2ml)e) 


The population build up is obtained by integrating 

r°° dfi j - 

n{eo,y)= de'—{e',y) ; sq = y"^ + mlr^ = yy/1 + uq . 


(B.3) 


where the initial population of Uh has been neglected and the value of Eq is set by the temper¬ 
ature at which pions appear in thermal equilibrium, assumed almost immediately after the 
hadronization transition. It is assumed that pion thermalization happens instantaneous ly at 
the QCD phase transition and we set Tq ~ 1 which is justified by the lattice results of 108 


which suggest a continuous transition that allows for thermalization on strong interaction 
time scales. 

The arguments of the exponentials inside of the logarithms can be shown to be < 0 so 
that the logarithms can be expanded consistently in a power series. Upon expanding the 
logarithms and integrating the rate, the population becomes 




6/J 


(B.4) 


where the expressions Xih{k), Jih{k) are given by 


Iih{k,y) = 
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(B.6) 
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+ (-1)^-^ / du- 

JuQ yUyl -|- u 
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■ I 1 I 


k yy/T+ 
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These expressions can be written in terms of the incomplete modified Bessel functions ini 
but this proves to be an unilluminating exercise. However, it is possible that slightly simpler 
expressions can be found under appropriate approximations such as the assumption that the 
neutrinos are produced either ultra-relativistically or non-relativistically over the duration 
of production. The dimensionless variable u = which was introduced in Eq IB. 11 


56 









































is the obvious quantity which controls how much of an ultra-relativistic or non-relativistic 
dispersion relation is obeyed by the heavy neutrino; specihcally u <C 1 implies a light, 
ultra-relativistic species while u ^ 1 implies a very heavy, non-relativistic species. 

From hgs 1^1171151 it can be seen that y, t take the values 1 < r < 10 and 0 < ?/ < 3 from the 
beginning of production until decoupling which implies that 0 < y/r < 3 for any neutrino 
produced. The production of purely ultra-relativistic neutrinos implies that u -C 1 or that 
rriy/m.,, <^y/T for all values of r, ?/ throughout production until freezeout. This assumption 
will hold for values of mvlmt^ -C 1 or for very large values of y. For ultra-relativistic 
production, IB.5IIB.6I simplify to 


Uk(k,y) - 

6/^ 
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1 y 1 + 


rrih) 
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2rn 


kA 


■ exp 


Ih 
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exp 
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—TV 


2m 
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(B.7) 


) \kAihy J \kA 




where r{u,z) is the incomplete gamma function. Further assuming Mh/M^^ -C 1 leads to 
another simplihcation 

^ (S|^) ^ 

With these approximations, the form of the distribution for production of purely ultra- 
relativistic neutrinos becomes 


where 


nih{ro,y) 


A 

^Ih 

'1+ 

exp 

-ky/Aih) 

UR 1/2 ^ 

^ k=l 

i + ey 
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Jk{'ro,y) 


(B.9) 
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F 1/2, 


kAihTl\ 

4// J 


(B.IO) 

This expression is equivalent to the distribution which was obtained in and explicitly 
depends on the time at which production begins, tq. In ^ it is shown that this approximate 
distribution is in excellent agreement with the distribution obtained from the exact numerical 
solution provided that Mh < IMeV. 

From hgs IbllTIIHl it is clear that for massive neutrinos with Mh > lOMeld there is a 


very differently shaped distribution function at small values of y co mpared with that seen 
from heavy neutrinos with masses below an MeV (as studied in [4^). This switch in shape 
between the distributions can be understood in terms of the production of non-relativistic 
neutrinos such that u 3> 1. The assumption that u 3> 1 implies that y/r for 

all values of r, i/ and, since < 1, this approximation is clearly only valid for a range 

of values for i/ -C 1 which explains the plateau in the distribution function for the heavier 
species. Under this assumption, the expressions IB.5IIB.6I become 
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where Ei{x) is the exponential integral. In the ultra-relativistic limit, exploiting 
was used in order to ignore several terms up to leading order in this allowed the 

sum, from lB.5lB.6l to be evaluated analytically. With heavier sterile neutrinos, such an 
approximation is unavailable and we are forced to retain the expression in the form of IB.Ill 
if we desire any nontrivial momentum dependence. 
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